On variational iterative methods for semilinear problems

TL;DR

This paper introduces an iterative approach for semilinear problems that simplifies them into linear systems solvable by fast Poisson solvers, demonstrating efficiency and accuracy through theoretical analysis and experiments.

Contribution

The paper proposes a novel iterative method that reduces semilinear problems to linear systems, enabling faster and accurate solutions compared to existing methods.

Findings

Method reduces computational cost significantly.

Achieves high accuracy in numerical approximations.

Supported by theoretical analysis and experimental validation.

Abstract

This paper presents an iterative method suitable for inverting semilinear problems which are important kernels in many numerical applications. The primary idea is to employ a parametrization that is able to reduce semilinear problems into linear systems which are solvable using fast Poisson solvers. Theoretical justifications are provided and supported by several experiments. Numerical results show that the method is not only computationally less expensive, but also yields accurate approximations.

Tables

Table 1. Table 1 : Relative error of series approximation for various M 𝑀 M and number of elements.

		Number of elements
	$2^3\times 2^3$	$2^4\times 2^4$	$2^5\times 2^5$	$2^6\times 2^6$	$2^7\times 2^7$
$M=1$	$3.80095\left(-01\right)$	$9.56161\left(-02\right)$	$2.99741\left(-02\right)$	$1.41233\left(-02\right)$	$1.02027\left(-02\right)$
$M=3$	$3.88708\left(-01\right)$	$9.61027\left(-02\right)$	$2.98165\left(-02\right)$	$1.38540\left(-02\right)$	$9.90701\left(-03\right)$
$M=5$	$3.89504\left(-01\right)$	$9.61117\left(-02\right)$	$2.98152\left(-02\right)$	$1.38520\left(-02\right)$	$9.90491\left(-03\right)$
$M=7$	$3.89612\left(-01\right)$	$9.61119\left(-02\right)$	$2.98152\left(-02\right)$	$1.38520\left(-02\right)$	$9.90491\left(-03\right)$
Picard	$3.89633\left(-01\right)$	$9.61119\left(-02\right)$	$2.98152\left(-02\right)$	$1.38520\left(-02\right)$	$9.90491\left(-03\right)$

Table 2. Table 2 : Relative error of series approximation for various M , N 𝑀 𝑁 M,N and number of elements.

$N=2$		Number of elements
	$2^3\times 2^3$	$2^4\times 2^4$	$2^5\times 2^5$	$2^6\times 2^6$	$2^7\times 2^7$
$M=1$	$1.77952\left(-01\right)$	$1.75391\left(-01\right)$	$1.72629\left(-01\right)$	$1.72591\left(-01\right)$	$1.72622\left(-01\right)$
$M=3$	$8.85117\left(-02\right)$	$3.72049\left(-02\right)$	$1.08392\left(-02\right)$	$5.66948\left(-03\right)$	$5.07559\left(-03\right)$
$M=5$	$8.90338\left(-02\right)$	$3.66827\left(-02\right)$	$9.31148\left(-03\right)$	$2.35028\left(-03\right)$	$6.04344\left(-04\right)$
$M=7$	$8.90409\left(-02\right)$	$3.66799\left(-02\right)$	$9.30586\left(-03\right)$	$2.34311\left(-03\right)$	$5.93845\left(-04\right)$
Picard	$8.90290\left(-02\right)$	$3.66693\left(-02\right)$	$9.29572\left(-03\right)$	$2.33310\left(-03\right)$	$5.83828\left(-04\right)$
$N=4$
$M=1$	$1.77952\left(-01\right)$	$1.75392\left(-01\right)$	$1.72632\left(-01\right)$	$1.72594\left(-01\right)$	$1.72624\left(-01\right)$
$M=3$	$8.84999\left(-02\right)$	$3.71949\left(-02\right)$	$1.08316\left(-02\right)$	$5.66720\left(-03\right)$	$5.07644\left(-03\right)$
$M=5$	$8.90218\left(-02\right)$	$3.66722\left(-02\right)$	$9.30139\left(-03\right)$	$2.34038\left(-03\right)$	$5.94682\left(-04\right)$
$M=7$	$8.90289\left(-02\right)$	$3.66693\left(-02\right)$	$9.29571\left(-03\right)$	$2.33314\left(-03\right)$	$5.83835\left(-04\right)$
Picard	$8.90290\left(-02\right)$	$3.66693\left(-02\right)$	$9.29572\left(-03\right)$	$2.33310\left(-03\right)$	$5.83828\left(-04\right)$

Table 3. Table 3 : Relative error of series approximation for various M , N 𝑀 𝑁 M,N and number of elements.

$N=2$		Number of elements
	$2^3\times 2^3$	$2^4\times 2^4$	$2^5\times 2^5$	$2^6\times 2^6$	$2^7\times 2^7$
$M=1$	$3.67061\left(-01\right)$	$3.25203\left(-01\right)$	$3.18984\left(-01\right)$	$3.17578\left(-01\right)$	$3.17235\left(-01\right)$
$M=3$	$1.02834\left(-02\right)$	$3.62913\left(-02\right)$	$2.20172\left(-02\right)$	$1.93606\left(-02\right)$	$1.88047\left(-02\right)$
$M=5$	$9.31325\left(-02\right)$	$2.41869\left(-02\right)$	$6.25023\left(-03\right)$	$1.74825\left(-03\right)$	$6.98383\left(-04\right)$
$M=7$	$9.29634\left(-02\right)$	$2.40161\left(-02\right)$	$6.06422\left(-03\right)$	$1.52422\left(-03\right)$	$3.86019\left(-04\right)$
Picard	$9.29576\left(-02\right)$	$2.40101\left(-02\right)$	$6.05810\left(-03\right)$	$1.51804\left(-03\right)$	$3.79729\left(-04\right)$
$N=4$
$M=1$	$3.67060\left(-01\right)$	$3.25202\left(-01\right)$	$3.18984\left(-01\right)$	$3.17578\left(-01\right)$	$3.17235\left(-01\right)$
$M=3$	$1.02830\left(-02\right)$	$3.62881\left(-02\right)$	$2.20151\left(-02\right)$	$1.93592\left(-02\right)$	$1.88035\left(-02\right)$
$M=5$	$9.31286\left(-02\right)$	$2.41828\left(-02\right)$	$6.24614\left(-03\right)$	$1.74434\left(-03\right)$	$6.95296\left(-04\right)$
$M=7$	$9.29585\left(-02\right)$	$2.40103\left(-02\right)$	$6.06010\left(-03\right)$	$1.52007\left(-03\right)$	$3.80026\left(-04\right)$
Picard	$9.29576\left(-02\right)$	$2.40101\left(-02\right)$	$6.05810\left(-03\right)$	$1.51804\left(-03\right)$	$3.79729\left(-04\right)$

Equations

$$-\nabla \cdot (a\nabla \psi )+r(\text{x},\psi )=f(\text{x}),\text{x}\in \Omega ,\psi (\text{x})=0,\text{x}\in \partial \Omega ,$$

$$H^{1}_{0}(\Omega)=\big{\{}w:w\in L^{2}(\Omega),~{}\nabla w\in L^{2}(\Omega)~{}\text{ and }~{}w|_{\partial\Omega}=0\big{\}}.$$

$$\begin{array}{rl}& \text{Find u∈H01(Ω) satisfying}\\ & {\int }_{\Omega }a\nabla u\cdot \nabla w\text{d}A+{\int }_{\Omega }r(\cdot ,u)w\text{d}A={\int }_{\Omega }fw\text{d}A,\forall w\in H_0^1(\Omega ).\end{array}$$

$$\mathcal{V}_{h}=\big{\{}w_{h}\in C(\overline{\Omega}):w_{h}|_{\tau}\text{ is linear for all }\tau\in\mathcal{T}_{h},\text{ and }w_{h}|_{\partial\Omega}=0\big{\}}.$$

$$\begin{array}{rl}& \text{Find uh∈Vh satisfying}\\ & {\int }_{\Omega }a\nabla u_h\cdot \nabla w_h\text{d}A+{\int }_{\Omega }r(\cdot ,u_h)w_h\text{d}A={\int }_{\Omega }f{w}_h\text{d}A,\forall w_h\in {\mathcal{V}}_h.\end{array}$$

$$\mathcal{V}_{h}=\text{span}\big{\{}\phi_{z}\big{\}}_{z\in Z}$$

$$\mathcal{R}(\mathbf{\alpha })=\mathcal{A}\mathbf{\alpha }+\mathcal{B}(\mathbf{\alpha })+\mathbf{\beta }=0,$$

$${\mathcal{A}}_{\zeta z}={\int }_{\Omega }a\nabla {\varphi }_z\cdot \nabla {\varphi }_{\zeta }\text{d}A,{\beta }_{\zeta }={\int }_{\Omega }f{\varphi }_{\zeta }\text{d}A\text{and}{\mathcal{B}}_{\zeta }(\mathbf{\alpha })={\int }_{\Omega }r(\sum {\alpha }_z{\varphi }_z){\varphi }_{\zeta }\text{d}A.$$

$$u_h=\sum _{z\in Z}{\alpha }_z{\varphi }_z.$$

$$\psi =\sum _{\zeta \le M}{\psi }_{\zeta }\text{and}r=\sum _{\zeta \le M}{\mathcal{P}}_{\zeta }$$

$$\mathcal{P}_{\zeta}=\dfrac{1}{(\zeta-1)!}\dfrac{\textrm{d}^{\zeta-1}}{\textrm{d}\lambda^{\zeta-1}}\Big{[}r\Big{(}\sum_{j\geq 0}\lambda^{j}\psi_{j}\Big{)}\Big{]}_{\lambda=0}.$$

$$\sum _{\zeta \le M}\left({\int }_{\Omega }a\nabla {\psi }_{\zeta }\cdot \nabla w\text{d}A+{\int }_{\Omega }{\mathcal{P}}_{\zeta }w\text{d}A\right)={\int }_{\Omega }fw\text{d}A.$$

$$\begin{array}{rl}& \text{Find ψζ∈H01(Ω) satisfying}\\ & {\int }_{\Omega }a\nabla {\psi }_{\zeta }\cdot \nabla w\text{d}A=({\delta }_{\zeta 0}-1){\int }_{\Omega }{\mathcal{P}}_{\zeta }w\text{d}A+{\delta }_{\zeta 0}{\int }_{\Omega }fw\text{d}A,\forall w\in H_0^1(\Omega ).\end{array}$$

$${\psi }_M=\sum _{\zeta \le M}{\psi }_{\zeta }.$$

$$\Vert \nabla {\psi }_{\zeta }{\Vert }_2\le {\delta }_{\zeta 0}{C}_{\zeta }\Vert f{\Vert }_2+(1-{\delta }_{\zeta 0}){\tilde{C}}_{\zeta }\Vert {\mathcal{P}}_{\zeta }{\Vert }_2,\zeta =0,1,2,\cdots ,M,$$

$$\Vert \nabla {\psi }_M{\Vert }_2\le C_{\zeta }\Vert f{\Vert }_2+\sum _{\zeta \le M}{\tilde{C}}_{\zeta }\Vert {\mathcal{P}}_{\zeta }{\Vert }_2.$$

$${\int }_{\Omega }a\nabla {\psi }_M\cdot \nabla w\text{d}A+\sum _{\zeta \le M}{\int }_{\Omega }{\mathcal{P}}_{\zeta }w\text{d}A={\int }_{\Omega }fw\text{d}A.$$

$$\begin{array}{rl}{\int }_{\Omega }a\nabla e_M\cdot \nabla w\text{d}A& ={\int }_{\Omega }r(\psi )w\text{d}A-\sum _{\zeta \le M}{\int }_{\Omega }{\mathcal{P}}_{\zeta }w\text{d}A=\sum _{\zeta >M}{\int }_{\Omega }{\mathcal{P}}_{\zeta }w\text{d}A.\end{array}$$

$$\begin{split}\|\nabla e_{M}\|^{2}_{2}&\leq\Bigg{(}\sum_{\zeta>M}\|\mathcal{P}_{\zeta}\|_{2}\Bigg{)}\|e_{M}\|_{2}\leq\Bigg{(}\sum_{\zeta>M}\|\mathcal{P}_{\zeta}\|_{2}\Bigg{)}\|\nabla e_{M}\|_{2}.\end{split}$$

$$\Vert \nabla e_M{\Vert }_2\le \sum _{\zeta >M}\Vert {\mathcal{P}}_{\zeta }{\Vert }_2.$$

$$\begin{array}{rl}{\int }_{\Omega }a\nabla {\psi }_M\cdot \nabla w\text{d}A& ={\int }_{\Omega }fw\text{d}A-\sum _{\zeta \le M}{\int }_{\Omega }{\mathcal{P}}_{\zeta }w\text{d}A\\ & ={\int }_{\Omega }fw\text{d}A-{\int }_{\Omega }{\psi }_{M-1}w\text{d}A.\end{array}$$

$${\int }_{\Omega }a\nabla {\psi }^{(M)}\cdot \nabla w\text{d}A={\int }_{\Omega }fw\text{d}A-{\int }_{\Omega }{\psi }^{(M-1)}w\text{d}A.$$

$$\begin{array}{rl}{\int }_{\Omega }a\nabla ({\psi }_M-{\psi }^{(M)})\cdot \nabla w\text{d}A& ={\int }_{\Omega }\left({\psi }^{(M-1)}-{\psi }_{M-1}\right)w\text{d}A.\end{array}$$

$$\mathcal{Y}=\log a,$$

$$a=\sum _{\xi \le N}{\mathcal{Y}}^{\xi }/\xi !.$$

$$\sum _{\xi \le N}\frac{1}{\xi !}{\int }_{\Omega }{\mathcal{Y}}^{\xi }\nabla {\psi }_{\zeta \xi }\cdot \nabla w\text{d}A=({\delta }_{\zeta 0}-1){\int }_{\Omega }{\mathcal{P}}_{\zeta }w\text{d}A+{\delta }_{\zeta 0}{\int }_{\Omega }fw\text{d}A.$$

$$\begin{array}{rl}& {\int }_{\Omega }\nabla {\psi }_{\zeta 0}\cdot \nabla w\text{d}A=({\delta }_{\zeta 0}-1){\int }_{\Omega }{\mathcal{P}}_{\zeta }w\text{d}A+{\delta }_{\zeta 0}{\int }_{\Omega }fw\text{d}A,\\ & {\int }_{\Omega }\nabla {\psi }_{\zeta \xi }\cdot \nabla w\text{d}A+\sum _{j=1}^{\xi }\frac{1}{j!}{\int }_{\Omega }{\mathcal{Y}}^j\nabla {\psi }_{\zeta ,\xi -j}\cdot \nabla w\text{d}A=0,\xi \ge 1.\end{array}$$

$${\psi }_{MN}=\sum _{\zeta \le M}\sum _{\xi \le N}{\psi }_{\zeta \xi }.$$

$$\Vert \nabla {\psi }_{\zeta 0}{\Vert }_2\le C_{\zeta 0}\left({\delta }_{\zeta 0}\Vert f{\Vert }_2+\Vert {\mathcal{P}}_{\zeta }{\Vert }_2\right)$$

$$\begin{array}{rl}\Vert \nabla {\psi }_{\zeta 1}{\Vert }_2& \le \Vert \mathcal{Y}{\Vert }_2\Vert \nabla {\psi }_{\zeta 0}{\Vert }_2\\ & \le C_{\zeta 1}\Vert \mathcal{Y}{\Vert }_2\left({\delta }_{\zeta 0}\Vert f{\Vert }_2+\Vert {\mathcal{P}}_{\zeta }{\Vert }_2\right).\end{array}$$

Full text

Keywords: Finite element method, Adomian polynomials, Picard iteration

1 Introduction

The object of discussion in this article is an iterative procedure for inverting the semilinear problem

[TABLE]

where $\Omega$ is a bounded region with boundary $\partial\Omega$. The functions $f$ and $g$ are known continuous data and $r:\Omega\times\mathbb{R}\to\mathbb{R}$ is continuously differentiable. In a variety of context, (1.1) is of great importance to mathematicians, physicists, biologists, engineers and other practitioners for modeling natural and human-induced phenomena. Models arising in this form describe steady-state characterization of evolution processes and have been the subject of extensive research in recent years. Some of the most commonly encountered classical mathematical models are prototypes of the differential equation in (1.1). They include the well-known Helmholtz equation, Bratu, Poisson-Boltzmann, the stationary Klein-Gordon, Allen-Cahn equation and Schrödinger equations. Precisely, applications arise in, but not limited to the study of electrostatic properties [5] and the study of acoustic and electromagnetic scattering [8, 16, 22]. In the sphere of fluid dynamics, some schemes (e.g. the Implicit Pressure, Explicit Saturation formulation of water-oil displacement [10]) are governed by variants of (1.1).

Indeed, simpler variations of (1.1) are separable in some coordinate systems, which allow solutions to be written in closed form. However, the fundamental models governing most natural phenomena involving motion, diffusion, reaction, etc. occur nonlinearly. Despite their pivotal contribution to the evolution of mathematical theories and experimentation, analytical methods have served minimal roles as only a handful of realistic models have known solutions. For several decades now, scientific computation has provided critical information about a wide class of problems that would otherwise be unavailable via analytical means. From a computational standpoint, solving (1.1) requires the simultaneous application of a numerical scheme with some sort of iterative procedure to an associated discrete algebraic system. As it is known, classical central and upwind schemes for the model problem (1.1) do not always produce satisfactory approximations. This is due to the inherent complex dynamics of $r$ leading to numerical challenges such as the appearance of singular effects [2, 19]. Even though these numerical tools are very successful, resolving instabilities for stable and accurate approximations often require a lot of computational time and large amounts of memory.

Ideally, analytical solutions to differential equations are desired. While this is possible when some components of the equation are trivialized, realistic models rarely admit closed form solutions. In recent years, advances in differential and integro-differential equations saw the light of novel analytical techniques [1, 3, 11, 23, 24]. The heart of these methods relies on perturbation arguments and the conjecture that solutions to nonlinear equations can be expressed as infinite series. This leads to a set of explicit recursive relations governing the modes of the series. Methods of this sort are robust and efficient for approximating solutions to IVPs, but have been modified to adapt to boundary and initial boundary-value problems. Certainly, the principal theme that has commanded most discussions about the aforementioned techniques is the rapid convergence of the series solutions. Despite this conception, important questions about the convergence of the method are still open because related studies are mainly problem-specific and do not provide a unifying ground to generalize results.

Far from a computational perspective is the issue of existence, solvability, uniqueness, regularity and number of solutions to (1.1). Nonlinear differential equations, in general, do not yield easily to rigorous and consolidated mathematical analysis. Successful investigations into the existence and uniqueness of solutions have relied on rigid assumption on $a,r$ and $f$. In the notes [17], Naito investigated the existence of entire solutions by enforcing separability of $r$ in the spatial coordinates and the unknown. In particular, the form $r=a(\textbf{x})g(w)$ was assumed, where $a$ is locally Hölder on $\mathbb{R}^{n}$ and $f$ is a locally Lipschitz continuous function on $\mathbb{R}^{+}$ which is positive and nondecreasing for $w>0$. In the said article, the author established approximate range of the limit values of positive solutions. Naito’s results were later improved by Usami [21], where the general form $r=g(\textbf{x},w,\nabla w)$ was considered. A more general case was studied by Boccardo and co. [6]. Their study proved the existence and regularity of (1.1) with a non-constant conductivity and $r=f(\textbf{x})w^{-\gamma}$, where $\gamma>0$. The article concluded that a locally and strictly positive distributional solution to (1.1) exists if $f\in L^{1}(\Omega)$. A related discussion that runs across the aforementioned studies (also see [2, 4, 14]) is the existence of multiple nontrivial solutions even under very restrictive assumptions. Sometimes, infinitely many solutions exist [13, 18]. As much as extensive study of special classes of nonlinear differential equations exist today, there is no unified theory of nonlinear equations [20].

Inspired by the monographs [3, 15], this study presents an iterative procedure for approximating the solution, if it exists, to the variational equation associated with (1.1). The method is developed by employing a decomposition that is able to reduce the problem into a collection of sub-equations governed by the Laplacian. For compatible boundary conditions and domains, the method allows Fast Poisson Solvers to be used. The method extends to other boundary types trivially. Prescribing pure Dirichlet or Neumann boundary conditions is a matter of convenience and allows Fourier transforms to be used.

2 Variational formulation

Consider a bounded Lipschitz domain $\Omega\subset\mathbb{R}^{2}$ with polygonal boundary $\partial\Omega$. Let $a:\Omega\to\mathbb{R}^{+}$ with $0<a_{\min}\leq a(\textbf{x})\leq a_{\max}$ for all $\textbf{x}\in\Omega$ and define the usual Hilbert space

[TABLE]

Under the assumption that a solution (not necessarily unique) exists, given $f\in L^{2}(\Omega)$, the variational formulation associated with (1.1) reads:

[TABLE]

In the sequel, the explicit dependence of $r$ on x is ignored and we write $r(\textbf{x},\psi)=r(\psi)$ for brevity. An important step towards the construction of the finite element solution consists of approximating the continuous problem in (2.3) by seeking the solution on a finite dimensional subspace of $H^{1}_{0}(\Omega)$. Mainly, $H^{1}_{0}(\Omega)$ is replaced by the set of piecewise continuous polynomials with compact support.

2.1 Galerkin approximation

Consider a triangulation $\mathcal{T}_{h}$ of $\Omega$ such that ${\overline{\Omega}}=\cup_{\tau\in{\mathcal{T}}_{h}}\tau$ with $h=\max_{\tau\in{\mathcal{T}}_{h}}\{h_{\tau}\}$ and $h_{\tau}=\sqrt{|\tau|}$. The finite element space associated with ${\mathcal{T}}_{h}$ is defined as

[TABLE]

Then, the associated Galerkin finite element formulation is:

[TABLE]

Let $Z$ be the set of degrees of freedom associated with $\mathcal{V}_{h}$ having $|Z|=\widetilde{N}$. Then,

[TABLE]

where $\phi_{z}$ is the finite element basis function tagged with the degree of freedom $z$. By setting $u_{h}=\sum_{z\in Z}\alpha_{z}\phi_{z}$, the variational problem in (2.5) reduces to the following discrete algebraic equation governing ${\boldsymbol{\alpha}}\in\mathbb{R}^{\widetilde{N}}$:

[TABLE]

where $\mathcal{A}\in M_{\widetilde{N}\times\widetilde{N}}(\mathbb{R}),~{}{\boldsymbol{\beta}},\in\mathbb{R}^{\widetilde{N}}$ and $\mathcal{B}({\boldsymbol{\alpha}})$ is a matrix function given by

[TABLE]

respectively. The Galerkin solution that approximates $\psi$ is written

[TABLE]

As $\mathcal{R}$ is nonlinear, direct inversion of (2.7) requires the use of some iterative scheme. Standard practice involves the application of quadrature rules alongside a linearization procedure. With appropriate regularity conditions, Newton methods guarantee higher convergence and are generally preferred. However, the laborious and unavoidable calculation of the associated Jacobian matrix makes Newton methods unappealing. Beside fixed point iterations, analytical methods have also been utilized to obtain qualitative features of solutions to nonlinear problems. Discussed in the next section is an iterative method which provides a general machinery for inverting (2.3).

3 An iterative approach

We furnish the weak form in (2.3) with a parametrization which transforms the variational problem into a sequential inversion of linear sub-problems. Suppose that $\psi$ and $r$ admit the decompositions

[TABLE]

where $\psi_{\zeta},~{}0\leq\zeta\leq M$, are to be determined. Here, construction of the set $\{\mathcal{P}_{\zeta}\}_{\zeta\leq M}$ is not well-defined, and thus, gives the reader the freedom to explore several possibilities. The construction adopted in this study engages sets of Adomian polynomials [11, 24] defined as

[TABLE]

This way, $\mathcal{P}_{\zeta}$ depends continuously on $\psi_{j},~{}0\leq j\leq\zeta-1$ only. Upon substituting (3.10) into (2.3), the variational equation can be written in the form

[TABLE]

The ultimate goal is to invert a set of linear sub-problems. Thus, the iterative scheme is developed by constructing the sequence $\{\psi_{\zeta}\}_{\zeta\leq M}$ where the terms in (3.12) are associated such that $\psi_{\zeta}$ solves

[TABLE]

where $\delta$ is the Kronecker delta. By construction, the term $\mathcal{P}_{\zeta}$ in the equation for $\psi_{\zeta}$ is known. Moreover, $\psi_{0}$ absorbs the boundary conditions of the original problem. All higher order terms serve as corrections to the approximation $\psi=\psi_{0}$ and honor the same boundary structure but homogeneous conditions. The practical solution is given by the $(M+1)$-term partial sum

[TABLE]

By requiring that $\mathcal{P}_{\zeta}\in L^{2}(\Omega)$ for all $0\leq\zeta\leq M$, Lax-Milgram supplies the existence and uniqueness of each mode, and hence the partial solution. Because the decomposition is performed on the variational equation, the approach is applicable to other boundary types with equal success. Choosing pure Dirichlet or Neumann condition are means to facilitate the use of Fast Fourier Transforms.

Remark 3.1**.**

It is noteworthy that, as it is, this formulation may yield the trivial partial solution. Observe that, for a homogeneous equation, the zeroth mode solves $\int_{\Omega}a\nabla\psi_{\zeta}\cdot\nabla w\,\text{d}A=0$ which yields the $\psi_{0}=0$ upon enforcing homogeneous boundary conditions. Largely, the zeroth mode inherits the most important qualitative features of the solution. As $\psi_{0}$ deviates significantly, convergence is affected adversely. This can be avoided by introducing an auxiliary function $g\in L^{2}(\Omega)$ into the differential equation. For instance, one may rewrite the original equation $-\nabla\cdot(a\nabla\psi)+r(\textbf{x},\psi)=0$ as $-\nabla\cdot(a\nabla\psi)+\widetilde{r}(\textbf{x},\psi)=g(\textbf{x})$, where $\widetilde{r}(x,\psi)=r(\textbf{x},\psi)+g(\textbf{x})$. This leads to a new set of polynomials $\{\widetilde{\mathcal{P}}_{\zeta}\}_{\zeta\leq M}$, and a nontrivial solution zeroth mode available by solving $\int_{\Omega}\nabla\psi_{0}\cdot\nabla w\,\text{d}A=\int_{\Omega}gw\,\text{d}A$.

Contrary to fixed point methods, this approximation eliminates any form of linearization. The following results provide theoretical justifications which support stability of the modes and also establish convergence of (3.14).

Lemma 3.2**.**

For all $\zeta\leq M$, the partial solution in (3.14) is stable provided $\mathcal{P}_{\zeta}\in L^{2}(\Omega)$ for all  $0\leq\zeta\leq M$.

Proof.

Stability of each mode is immediate since the weak forms in (3.13) possess unique solutions in $H^{1}_{0}(\Omega)$ with

[TABLE]

where $\|\cdot\|_{2}=\left(\int_{\Omega}|\cdot|^{2}\,\text{d}A\right)^{1/2}$ and the constants $C_{\zeta},\widetilde{C}_{\zeta},~{}\zeta=0,1,2\cdots$ are independent of $\psi_{\zeta}$. It follows that

[TABLE]

This completes the proof. ∎

Lemma 3.3**.**

Let $\psi$ satisfy (2.3) and $\psi_{M}$ be the partial solution in (3.14). If $\mathcal{P}_{\zeta}\in L^{2}(\Omega)$ for all $0\leq\zeta\leq M$, then $\|\nabla(\psi-\psi_{M})\|_{2}\to 0$ as $M\to\infty$.

Proof.

By assumption, $\psi$ solves (2.3) and $\psi_{M}$ satisfies

[TABLE]

It follows that the error $e_{M}=\psi-\psi_{M}$ satisfies

[TABLE]

From Cauchy-Schwarz, Poincaré and triangle inequalities, it holds that

[TABLE]

Consequently,

[TABLE]

The theorem follows as $\sum_{\zeta>M}\mathcal{P}_{\zeta}\to 0$ for sufficiently large $M$. ∎

3.1 Comparison with Picard iteration

The iterative method discussed above shares some similarities with Picard Iteration, where each iteration solves a linear problem. For the form in (2.3), the fixed point iteration is:

In fact, Lemma 3.4 shows that if $r$ is linear and the set $\{\mathcal{P}_{\zeta}\}_{\zeta\geq 0}$ is constructed according to (3.11), the two schemes are identical.

Lemma 3.4**.**

Let $\psi_{M}$ and $\psi^{(M)}$ be the $M^{\text{th}}$ iterate of the procedure above and Picard iteration respectively. Assume that $r$ is linear and $\{\mathcal{P}_{\zeta}\}_{\zeta\geq 0}$ is constructed according to (3.11). If $\psi_{0}=\psi^{(0)}$, then $\psi_{M}=\psi^{(M)}$ holds for all $M\geq 1$.

Proof.

The argument is inductive. Since $r$ is linear, we must have $\mathcal{P}_{\zeta}=r(\psi_{\zeta})=\psi_{\zeta}$. Thus,

[TABLE]

From algorithm 1, $\psi^{(M)}$ also satisfies

[TABLE]

Suppose for the sake of induction that $\psi^{(M-1)}=\psi_{M-1}$ for $M\geq 2$. It follows from (3.17) and (3.16) that

[TABLE]

The last equation is zero by assumption. As this holds for arbitrary $w\in H^{1}_{0}(\Omega)$, it must be that $\psi^{(M)}-\psi_{M}$ is constant in $H^{1}_{0}(\Omega)$. We then deduce that $\psi^{(M)}-\psi_{M}=0$. ∎

Equation (3.15) shows that convergence of the partial solution is controlled by the rate of decay of $\mathcal{P}_{\zeta}$. Certainly, this procedure is attractive because, similar to Picard iteration, calculation of the Jacobian associated with (2.7) is replaced by the repeated inversion of the operator $\nabla\cdot(a\nabla\cdot)$. Since the computational cost of inverting a single instance of this operator on a highly resolved mesh can be enormous for traditional numerical methods, calculation of (3.14) is highly prohibitive for large $M$. More importantly, engagement of Fast Poisson Solvers is limited because symmetry of the resulting stiffness matrix coefficients is affected.

In what follows, we reduce (3.13) into variational sub-problems governed by $\widetilde{a}(\textbf{x})=1$ for all x. This allows us to employ discrete sine transforms to accelerate calculation of the modes. To accomplish this, consider the conductivity log characterization [9]

[TABLE]

which exists, and is defined since $a(\textbf{x})>0$ for all $\textbf{x}\in\Omega$. This suggests that for sufficiently large $N$, $a$ can be approximated using the truncated Taylor series

[TABLE]

Using this representation for a fixed $\zeta$, the variational equation in (3.13) becomes

[TABLE]

By collecting terms of the same order (index of $\psi_{\zeta}$ and the exponent of $\mathcal{Y}$ adds up to the same $n\in\mathbb{N}$), (3.21) transforms into the following recursive equations:

[TABLE]

The approximate solution is written as

[TABLE]

Clearly, any unfavorable effect of $a$ on the invertibility of the resulting stiffness matrix is avoided. A minimal penalty is realized in the calculation of the load vector.

3.2 Discrete algorithm

The procedure discussed in the preceding sections can be summarized as follows:

For all $\zeta$, we can verify that the variational problems in Algorithm 2 have unique solutions in $H^{1}_{0}(\Omega)$ satisfying the following stability estimate:

[TABLE]

provided $\mathcal{P}_{\zeta}\in L^{2}(\Omega)$. By construction, $\int_{\Omega}\nabla\psi_{\zeta 1}\cdot\nabla w\,\text{d}A+\int_{\Omega}\mathcal{Y}\nabla\psi_{\zeta 0}\cdot\nabla w\,\text{d}A=0$. Thus, as long as $\mathcal{Y}\in L^{2}(\Omega)$, it follows that

[TABLE]

By complete induction, we have

[TABLE]

Upon a successful Galerkin discretization, the discrete equations have the form $\mathcal{M}{\boldsymbol{z}}=\mathcal{F}$ where $\mathcal{M}$ is the discrete Laplacian. The collection can now be solved efficiently and fast using the fast Poisson solver discussed in the next section. As much as the summation in the equation for $\psi_{\zeta\xi}$ tells a strong dependence of higher-order terms on lower-order ones, Fast Poisson solvers comes in handy especially when a lot of terms are required to represent each mode accurately.

For Neumann problems, the Galerkin discretization is known to produce a linear system characterized by a one-dimensional kernel. The associated stiffness matrix is singular with fatal effects on iterative solvers, e.g., conjugate gradient. Often, the kernel is eliminated by specifying data at a particular node, or approximate the system using a minimization-based iterative solver [7]. While these techniques seemingly avert direct inversion a singular matrix, Bochev and co., in the referenced article, remarked that such choices are subtle and unsatisfactory because specifying data at a specific node is ambivalent and roundoff errors may hinder convergence of the iterative solver. In this context, Algorithm 2 provides the environment to exploit direct solvers.

In very limited cases, (1.1) is fast solvable on simple geometries. The most common numerical treatment of the finite element discretization involves a uniform mesh refinement across the entire computational domain. This poses numerical challenges as the finite element method (as well as finite difference and finite volume) performs poorly with vanishing element width. For the prescribed boundary conditions, this is no longer a limitation because we are able to use Fourier transforms to invert the sub-systems efficiently.

4 A Fast Poisson Solver (FPS)

Fast Poisson Solvers are simple but powerful machinery which harness Fourier transforms to produce efficient solvers for the Poisson equation with compatible boundary conditions. They require a low amount of arithmetic operations and can be executed effectively with optimal operation count of $\mathcal{O}(n^{2}\log n)$. The procedure that utilizes finite difference is discussed in [12]. Let $\mathcal{M}{\boldsymbol{z}}=\mathcal{F}$ be the algebraic system resulting from the Galerkin approximation of

[TABLE]

using the bilinear basis. Then $\mathcal{M}$ consists of block diagonal matrices identified as TST (Toeplitz, Symmetric and Tridiagonal). Denote by TD$(n;\alpha^{\star},\beta^{\star})$ the $n\times n$ tridiagonal matrix with $\alpha^{\star}$ along the principal diagonal and $\beta^{\star}$ along the sub- and super-diagonals. For the finite element discretization, let $\Delta x$ and $\Delta y$ be the element width with $n_{1}$ and $n_{2}$ elements in the vertical and lateral directions respectively. Then $\mathcal{M}=$ TD$(n;\mathcal{S},\mathcal{T})$, where $\mathcal{S}$ and $\mathcal{T}$ are themselves TST with sub- and super-diagonals

[TABLE]

Interestingly, TST matrices encode critical information which are fundamental to the implementation of FFT-based Poisson solvers. In fact, the linear system $\mathcal{M}{\boldsymbol{z}}=\mathcal{F}$ is solvable completely by using the eigenvalues and eigenfunctions of $\mathcal{S}$ and $\mathcal{T}$ only. First, the eigenvectors are given explicitly by the formula

[TABLE]

Moreover, the collection $\{\varphi_{jl}\}_{jl=1}^{n_{1}}$ is orthonormal with corresponding eigenvalues

[TABLE]

With the set $\{\varphi_{jl},\lambda_{j}^{\mathcal{S}},\lambda_{j}^{\mathcal{T}}\}_{jl=1}^{n_{1}}$ available, construction of the fast Poisson solver is straightforward. Let $\boldsymbol{z}=(\boldsymbol{z}_{1},\boldsymbol{z}_{2},\cdots,\boldsymbol{z}_{n_{2}})^{T}$ and $\mathcal{F}=(\mathcal{Q}_{1},\mathcal{Q}_{2},\cdots,\mathcal{Q}_{n_{2}})^{T}$ where $\boldsymbol{z}_{l}$ and $\mathcal{Q}_{l}$ are block vectors corresponding to the $l^{\text{th}}$ column of $\mathcal{M}$ and $\mathcal{F}$ respectively. Then $\mathcal{M}{\boldsymbol{z}}=\mathcal{F}$ is expressible as

[TABLE]

Denote by $\Phi$ the matrix whose $(jl)^{\text{th}}$ entry is given by (4.26). Then

[TABLE]

where $D_{\mathcal{I}}=\text{diag}[\lambda_{l}^{\mathcal{I}},\lambda_{l}^{\mathcal{I}},\cdots,\lambda_{l}^{\mathcal{I}}]$. Orthogonality of $\Phi$ implies that (4.28) can be reduced to

[TABLE]

where, $\Phi\boldsymbol{z}_{l}=\boldsymbol{\zeta}_{l}$ and $\Phi\mathcal{Q}_{l}=\boldsymbol{\xi}_{l}$. These operations successfully decoupled $\mathcal{M}{\boldsymbol{z}}=\mathcal{F}$ into a set of independent one-dimensional tridiagonal systems which can be solved using fast tridiagonal solvers.

5 Numerical Tests

We present three numerical tests to explore the effect of $M$ and $N$ on the accuracy of the series approximations. Even though the problems are manufactured and may not bear any physical meaning, they do provide good insights into the performance of the method. Because of the close similarity, we also compute reference approximations by direct simulation of $\mathcal{R}({\boldsymbol{\alpha}})={\boldsymbol{0}}$ using Picard iteration. The stopping criterion is $\|\mathcal{R}({\boldsymbol{\alpha}})\|_{\infty}\leq 10^{-10}$.

5.1 Test 1

: In this foremost experiment, we set $a(x_{1},x_{2})=1,~{}r(\psi)=\psi-\psi^{3}$ and $f$ is chosen appropriately so that (1.1) admits the solution $\psi(x_{1},x_{2})=\omega(x_{1})\omega(x_{2})$, where

[TABLE]

Since $a$ is constant, the index $\xi$ is redundant. Thus, we simply denote $\psi_{\zeta\xi}$ by $\psi_{\zeta}$. Figure 1 compares the numerical approximations of $\psi_{1},\psi_{4}$ and $\psi_{8}$ with the analytical solution.

5.2 Test 2

: Here, we set $a(x_{1},x_{2})=\cosh(x_{1}+x_{2}),~{}r(\psi)=\sin\psi$ and $f$ is constructed such that (1.1) admits a closed form solution

[TABLE]

5.3 Test 3

: We set $a(x_{1},x_{2})=1+x_{1}^{2}+x_{2}^{2},~{}r(\psi)=\exp(-\psi^{2})+\tan^{-1}\psi$ and seek the closed form solution $\psi(x_{1},x_{2})=\sin(3\pi x_{1})\sin(2\pi x_{2})$. Moreover, $f$ is chosen accordingly with $r(\psi)=\psi$.

6 Discussions and conclusion

The numerical results show an impressive agreement between the series approximation and the exact solutions. Since the lower order modes do not inherit sufficient spatial variability from $a$, their approximations are distinguishable from the analytical solutions. Once the variability is well resolved, it takes only a few terms to obtain an accurate approximation. Moreover, the contribution of larger values of $M$ and $N$ to the overall error is insignificant as the finite element error dominates. The errors recorded show that using an excessive number of terms to construct the partial solution is obviously an overkill because the gain is insignificant compared to the computational cost. Overall, the method provides accurate approximations and the variational context provides a flexible framework to study its convergence. For the benchmark problem, the method allows the inversion of nonlinear differential operators to be achieved through a series of inversions of the discrete Laplacian. This is appealing because the nontrivial task of calculating the associated Jacobian matrix is avoided entirely. For linear reactions, the method is shown to coincide with the Picard iteration method. Even though the equivalence has not been established for a general nonlinear reaction, the numerical results of the two approaches are largely consistent.