On iterative methods for some elliptic equations with nonlocal conditions

Olga Štikonienėa,b,1, Mifodijus Sapagovasb,2, Regimantas Čiupaila Faculty of Mathematics and Informatics, Vilnius University Naugarduko str. 24, LT-03225 Vilnius, Lithuania olga.stikoniene@mif.vu.lt Institute of Mathematics and Informatics, Vilnius University Akademijos str. 4, LT-08663 Vilnius, Lithuania mifodijus.sapagovas@mii.vu.lt Vilnius Gediminas Technical University Saulėtekio ave. 11, LT-10223 Vilnius, Lithuania regimantas.ciupaila@vgtu.lt


Introduction and problem statement
Boundary value problems for various differential equations with nonlocal boundary conditions were actively investigated during the last three decades.A lot of attention was given to numerical analysis of such problems.Research is motivated by both their interest to pure mathematics and new applications in physics, mechanics, biochemistry, ecology.The first articles dealing with parabolic equations with integral conditions were related to heat equation [1] and the theory of linear thermoelasticity [2].Many other applications of problems with nonlocal boundary conditions can be found in books [3,4] while numerical methods are reviewed in papers [5][6][7][8][9][10].
One of the first papers, in which there were the new type boundary problems with nonlocal conditions for the elliptic equations formulated and investigated, was article [11].The nonlocal conditions introduced in this paper are the so called Bitsadze-Samarskii nonlocal conditions.The finite difference methods for the solution of the elliptic equations with the various nonlocal conditions were investigated in many works.However, these methods were rarely investigated for the parabolic equations.
One of the first papers in this particular area was article [12].In that paper, the finite difference method for the linear elliptic equation with the Bitsadze-Samarskii nonlocal condition was investigated.Latter, in the papers [13][14][15][16][17][18][19], the issues of the existence, uniqueness and convergence of the solution of the difference problem for the elliptic equations with the various types of the nonlocal conditions, including the integral ones, were considered.The recent results and the bibliography may be found in the works [20][21][22][23][24][25][26].
In the paper [23], the solution of the elliptic equation with the nonlocal conditions is also investigated bringing this problem to the problem with the Dirichlet condition.
The idea to reduce the system of difference equations with nonlocal conditions to the two systems of the same type with the Dirichlet conditions was also used earlier for the solution of the elliptic equations with the nonlocal conditions of the special type [27,28].But in all quoted papers, the iterative methods for solutions of the systems of difference equations with nonlocal conditions were not considered.
Iterative methods for finite difference systems obtained from elliptic equations with nonlocal boundary conditions are insufficiently analysed.Matrices of these systems of equations have a lot of specific properties.These matrices are usually nonsymmetric, and their spectrum is relatively complicated and depends on parameters and functions of nonlocal conditions.
Relatively simple one dimensional eigenvalue problem was investigated in article [28].There were noticed that the matrix of the difference equation can have zero, negative, complex or multiple eigenvalues depending on the values of γ and ξ.In the case of 2D elliptic problem, the spectrum of matrices of finite difference schemes is even more complicated and affects the convergence of iterative methods.We investigate a second order linear elliptic equation with the Dirichlet and nonlocal conditions where D = (0, 1) × (0, 1), p(x, y) > 0, q(x, y) 0. The conditions for the weight function α(x) will be formulated later.
Main aim of the paper is to investigate various methods for the solution of the system of difference equations, corresponding to the differential problem (1)-( 3).The convergence of few simple iterative methods could be justified, requiring that function α(x) should satisfy the conditions, guaranteeing that the matrix of the system of difference equations is M-matrix.When p(x, y) = const, the convergence of more effective iterative methods could be proven by investigating the structure of the spectrum of the matrix of the system of difference equations.This matrix is nonsymmetric because of the nonlocal condition (2), so, the structure of its spectrum might be more complicated than for the Dirichlet conditions.
It is significant that, in the case of λ = 0, solution (for boundary value problem) is not unique; λ < 0 (or λ ∈ C, Reλ < 0) may lead to instability of the difference scheme for the parabolic problem or to divergence of the iteration method for the stationary problem.
The remaining part of this paper is organized as follows.In Section 2, we state the difference problem.In Section 3, we analyse when the matrix of difference problem with nonlocal conditions is a M-matrix.The structure of the spectrum of difference problem is investigated in Section 4. In Sections 5 and 6, the convergence of the Chebyshev method and the alternating direction method is considered.The results of numerical experiments are presented in Section 7. Section 8 contains some brief conclusions and comments.

Finite-difference scheme
In the domain D, we introduce the grids: In D, we consider the rectangular grids If ω is one of these grids, we define the space of grid functions F(ω).We introduce first and second order difference operators: The functions f , q are approximated by grid functions f ij and q ij on the grid ω h , (µ 1 ) j , (µ 2 ) j on the grid ω h y , (µ 3 ) i , (µ 4 ) i on the grid ω h x .
Problem (1)-(3) in D is replaced by a finite-difference scheme with nonlocal condition The truncation error of the scheme is O(h 2 ).

Problems with nonlocal conditions and M-matrices
Equations ( 4)-( 6) can be rewritten in a matrix form where A = {a kl } -square matrix, k, l = 1, n, n = N (N − 1), a kk > 0. Now we give a short summary of useful results from the theory of M-matrices.More details can be found in [29,30].The most important properties of M-matrices are the discrete maximum principle and the guaranteed convergence of various linear solvers.
Definition.The square matrix A = {a kl }, k, l = 1, n, is called M-matrix if a kl 0 when k = l and there exists the inverse matrix, whose all elements are nonnegative (A −1 0).
Lemma 1.If the elements of the matrix A = {a kl } have a property a kl 0 when k = l, then three sequent statements are equivalent: (i) A −1 exists and A −1 0; (ii) real parts of all the eigenvalues of the matrix A are positive: Re λ(A) > 0; (iii) matrix A can be expressed in the form where S 0, a > (S).
So, each of the properties along with the property a kl 0 when k = l could be treated as a definition of the M-matrix.Lemma 2. If A is M-matrix, then necessarily a kk > 0.
Lemma 3. If a kk > 0, a kl 0 when k = l and matrix A is strongly diagonally dominant, i.e.
Lemma 4. If a kk > 0, a kl 0 when k = l and matrix A is a weakly diagonally dominant and irreducible, then matrix A is M-matrix.
Remark.The diagonal dominance is not an necessary condition for matrix A to be M-matrix.

Lemma 5. If
A is an M-matrix and D is a diagonal matrix with the nonnegative elements, then A + D is an M-matrix.Lemma 6.If A is an M-matrix, then a matrix B, obtained from the matrix A by equating to zero any non-diagonal element, is also an M-matrix.
Suppose, A is an M-matrix and A has a representation A = M−N, where M −1 0 and N 0. In this case, A has a regular splitting A = M − N. One of these splittings is referred in Lemma 1(iii): A = aI − S, where S 0, a > (S).

Lemma 7. If A is an M-matrix and A has a regular splitting
Corollary.If A is an M-matrix and A has a regular splitting A = M − N, then, for the system Ax = f , the iterative method The following statement will be proven for the system of difference equations ( 4)-( 6), rewritten in the matrix form (7). Before that, the auxiliary statement on the trapezoid rule for the numerical integration will be formulated.
then, for sufficiently small values of h, the inequality is valid.
The statement of Lemma 8 follows from properties (8) of the function α(x) and the error estimation for trapezoid rule Theorem 1.If p(x, y) > 0, q(x, y) 0 and the conditions of Lemma 8 for the function α(x) are satisfied, then, for all sufficiently small h, the matrix of system ( 7) is an M-matrix.
Proof.The matrix A is weakly diagonally dominant.The off-diagonal entries of this matrix are nonpositive.Furthermore, the matrix of the system of difference equations (4) in rectangle is irreducible.So, according to Lemma 4, the matrix A of system ( 7) is an M-matrix.Now we can provide a few converge iterative methods for the solution of system ( 7), taking different regular splitting of the matrix A. With that aim, we will use the expression of the matrix where D, −L and −R are diagonal, strictly lower triangular and strictly upper triangular parts of A, respectively.Obviously, D > 0, L 0, R 0.
According to this expression, we define three different regular splittings In the case And the splitting with All these iterative methods converge under the conditions of Theorem 1.
We should admit that condition α(x) 0 in ( 8) is a necessary condition for the matrix A to be a M-matrix.In contrast, nondiagonal elements in equations ( 6) would be of different signs.
For one-dimensional parabolic equations, a bit different condition [31] could be formulated instead of conditions (8): For problem (1)-(3) and condition (10), we investigated that the conclusion that the matrix of system (7) would be an M-matrix does not follow.But this matrix is diagonally dominant under condition (10).This is enough for the reasoning of the methods of Jacobi, Gauss-Seidel and fixed point iterations.
As it follows from conditions ( 8) and (10), in the case of nonconstant α(x), max α(x) could be a sufficiently big number.But in the case when α(x) = const, it follows from (8) that Theorem 1 is true only under condition |α| < 1.
In the next section, we will demonstrate that, in the case of problem (1)-(3) when α(x) is a constant, the limitation α < 1 is not necessarily condition for the matrix A of system (7) to be an M-matrix.To this end, we consider the eigenvalue problem for the difference operator with the nonlocal condition.

Structure of the spectrum of the difference operator
Let us take the constant coefficients in problem ( 1)-(3), i.e. p(x, y) = 1, q(x, y) = 0, α(x) = γ 0. So, we consider the differential problem The corresponding difference problem could be written as follows: Let us express U 0j from equations (15): where Substituting the expression of U 0j into the equations ( 14), when i = 1, we get hα with the boundary conditions (16).For expression (17) to be always valid, it is enough to assume that h depending on γ would be sufficiently small, i.e.
Under that presumption, we get α 0. Let us change the sign on the both sides of equations ( 18), (19) and rewrite this system in a matrix form This system of the equations is composed analogically as system (7).The main difference is that the order of the matrix A 1 is (N − 1) 2 instead of N (N − 1) and the vector U is constructed from the values U ij taken from the internal points of the domain ω h .By solving this system, the values of U 0j , j = 1, N − 1, could be obtained according to formulas (17).

Lemma 9. The eigenvalue problem of the matrix
is an equivalent problem for the difference problem with the nonlocal conditions Proof.We express from equations (23).We get (21) substituting these expressions into equation ( 23), when i = 1.Now, let us consider problem (21).Let us define the values U 0j , j = 1, N − 1, by formulas (25).Then we could rewrite problem (21) using expressions (25) in the form of ( 22)- (24).The equivalence of the statement is proven.
Corollary.The number of the eigenvalues of the difference problem with the nonlocal conditions (22)-( 24) is (N − 1) 2 .
Now we should investigate, what are the conditions required for the matrix A 1 to be an M-matrix.The matrix A 1 as well as A has the property: a kl 0, when k = l, if γ 0 and h < 2/γ.To answer the question when would A 1 be an M-matrix, we should investigate the eigenvalue problem ( 22)- (24), or in other words, we should investigate the spectrum of the matrix A 1 .
Let us write the discrete solution of problem ( 22)-( 24) as So, we consider two auxiliary problems: (I) Problem with the nonlocal boundary conditions (II) Classical problem The eigenvalues of problem ( 22)-( 24) can be written as here are the eigenvalues for the second auxiliary problem ( 28)-( 29).
The eigenvalue problem ( 26)-( 27) has been considered in [32,33].We can formulate the main results of these papers in the following lemma.
From expressions (30) and Lemma 10 we derive the following statement.Now we can formulate and prove the theorem on the condition when the matrix A 1 is an M-matrix.
Theorem 2. There exists such a number γ 0 > 2 that, for all the numbers 0 γ γ 0 , the matrix A 1 is an M-matrix.
Proof.Suppose 0 γ < 2. Then all the eigenvalues η k are positive according to Lemma 10.The eigenvalues µ l are positive independent of γ.There exists only one eigenvalue η 1 = 0 when γ = 2. So, λ kl > 0 if 0 γ < 2. When γ > 2, just one and only one eigenvalue η 1 is negative The number β is the unique positive root of equation ( 31) that grows up as γ grows.Let us indicate the value of β 0 as the value of β, when µ 1 = |η 1 |, i.e.
Remark.For sufficiently small h, the value of β 0 described by formula (32) could be written approximately β 0 ≈ π.
Let us now we take a more general equation instead of equation ( 11) Let us consider the difference problem for the boundary value problem (34), ( 12), ( 13), analogous to the difference problem ( 14)- (16).We can write the system of these difference equations in the matrix form It is easy to see that matrix A 2 could be obtained from the matrix A 1 by adding q ij to the diagonal elements of the matrix A 1 .So, according to Lemma 5 and Theorem 2, the following statement is true.
Theorem 3.There exists such a number γ 0 > 2 (defined by formula (33)) that, for all values of 0 γ γ 0 , the matrix A 2 is an M-matrix.

Chebyshev method
In Sections 3 and 4, the systems of difference equations ( 7) and ( 21) were considered.Both the systems were contained with the specific feature, the matrix of the system was Mmatrix.This property allowed to construct the convergent iterative methods.We should note that the following condition was necessary to fulfill, the function α(x) in the nonlocal condition (2) and the parameter γ in analogous condition (12) should be nonnegative.
In this and the following sections, we will emphasize another property of the matrix of the system of difference equations, all the eigenvalues of the matrix of the system are positive, and the matrix itself is diagonalizable (the system of the eigenvectors of this matrix is complete or, in other words, the matrix is of the simple structure).Precisely, this property, along with certain conditions, is typical for many systems of difference equations, arising from the elliptic equation with various types of the nonlocal conditions.We admit that the nonlocal integral conditions derive that the matrix of the system of difference equations ( 21) is always nonsymmetric.
, where γ 0 is defined by formula (33), then all the eigenvalues of matrix A 1 of the system of difference equations (21) are positive, and the corresponding eigenvectors are linearly independent.
Proof.It follows from definition (30) of the eigenvalues of the matrix A 1 , Lemma 10 (when γ < 2) and Theorem 2 (when 0 γ < γ 0 ) that all the eigenvalues of the matrix A 1 are positive.Furthermore, the eigenvalue problems ( 26)-( 27) and ( 28)-( 29) are characterized by the feature that their systems of the eigenvectors are complete in the vector space R n−1 .Then it follows from the equality that the eigenvectors U k,l of the matrix A 1 are linearly independent (for more details, see [20]).We define iteration where The iterative method ( 36) is usually referred to as the Chebyshev method or the Richardson method, the set τ k , k = 1, m, is a Chebyshev set of iterative parameters.In the case of a symmetric positive defined matrix with the eigenvalues in [a, b], 0 < a < b, the Chebyshev iteration is optimal over all other polynomial based methods.The iterative method has the best convergence rate after m 1 iterations.
For the estimation of the convergence rate of the Chebyshev iterative method (36) in the case of a difference system with the nonlocal conditions, it is necessary to introduce norms for the vectors and the matrices.To this end, we define the matrix norm in a special way [6,34].By P we denote a matrix whose columns are linearly independent eigenvectors U k,l of the matrix A 1 (of the eigenvalue problem ( 22)-( 24)).We define a norm of an arbitrary square matrix B of order (N − 1) 2 by formula

For this definition, we have
A 1 * = (A 1 ).The norm of (N −1) 2 -dimensional vector U, compatible with this matrix norm, is defined as follows: where U = (U, U) 1/2 .Note that where (PP) −1 is a symmetric positive definite matrix.In this sense, the norm U * is the same as the energy norm considered in [35].
The following theorem is true.
Theorem 4. If −∞ < γ γ 0 , where γ 0 is defined by (33), then the Chebyshev iterative method (36) is convergent and the next estimate is true: where and U is the exact solution of system (21).
The theorem was proved under more general conditions in [34].As far as we are aware, the Chebyshev iterative method for the system of the difference equations with the nonlocal boundary conditions was considered for the first time in papers [34,36].The Chebyshev method for solution of the system of the linear equations with nonsymmetric matrix was investigated early in [37].

Peaceman-Rachford ADI method
The study of the structure of the spectrum of the difference operators with the nonlocal boundary conditions can be successfully applied in the investigation of convergence of alternating direction methods.
For each fixed value j = 1, N − 1, the system of equations ( 37) is solved with the nonlocal boundary condition and condition Analogously, for each fixed value i = 1, N − 1, system ( 38) is solved with the Dirichlet boundary conditions We rewrite the iterative method ( 37)- (38) with the corresponding boundary conditions in the matrix form.For this end, we introduce two (N − 1)-order matrices , where α = γ/(1 − hγ/2), and . Now we can definite two matrices Ã1 and Ã2 of the order (N − 1) 2 , using the Kronecker (tensor) product of matrices where I N −1 is an identity matrix of the order (N − 1).The iterative Peaceman-Rachford method ( 37)-( 41) can be written in matrix form: where Lemma 14.For the eigenvalues of the matrix S, the following equality is valid Lemma 15.If λ( Ã1 ) = a + ib and a > 0, τ k+1 > 0, then The validity of next theorem follows from Lemma 10 and Lemma 15.
We note now that, in the case γ < 2, all eigenvalues λ i ( Ã1 ) and λ j ( Ã2 ), i, j = 1, N − 1, are positive.But this limitation not necessary.It follows from Lemma 15 that the Peaceman-Rachford method converges then the difference operator with the nonlocal conditions has complex eigenvalues with positive real part.Such situation is possible for the differential equation (11) with the nonlocal condition (2) [39].
If there exists λ i ( Ã1 ) < 0, we give for λ i ( Ã1 ) < 0 under a condition that |λ i ( Ã1 )| is small enough.Such situation we observed in the numerical experiment.The same questions of the Peaceman-Rachford ADI method for the two-dimensional differential equation with the integral conditions were considered in [40][41][42][43].Various modifications of the alternating direction method and the locally one-dimensional methods were used in [44][45][46] for parabolic equations with the nonlocal conditions.

Numerical experiment
We consider a model problem . We set a reaction term q(x, y) = 1 and a weight function α(x) = 1.
A source function f (x, y), the initial and boundary conditions were prescribed to satisfy a given exact solution u(x, y) = 1 + e x+y .
We consider the uniform grids with the different mesh sizes h and analyze the convergence and the accuracy of the computed solution from the present ADI scheme.We compute the maximum norm of the error of the numerical solution with the respect to the exact solution, which is defined as We define a number p as p = ε 2h ε h , which theoretically must be approximately p ≈ 4.
Test problems were solved with the different values of the parameters γ and h.The numerical tests were carried out on the grids 2 k × 2 k , k = 2, . . ., 6.The computational results are reported in Tables 1 and 2.
We can clearly observe a second-order convergence in the maximum norm for all choices of γ 0 and γ 1 .We note that the Peaceman-Rachford method converges (Table 2) also for some cases λ j ( Ã1 ) < 0. The set of optimal iterative parameters of the ADI method has been chosen according to the monograph [47], where the symmetric matrices of an iterative process were used.

Conclusions and remarks
The iterative methods can be used for a numerical solution of the elliptic equations with the nonlocal conditions.Under some conditions on weight function α(x), the matrix of the system of difference equations is M-matrix.A linear system with an M-matrix can be efficiently solved by basic iterative methods, for example, by Gauss-Seidel method.When p(x, y) = const, the convergence of faster iterative methods was proven by investigating the structure of the spectrum of the matrix of the corresponding difference problem.This matrix is nonsymmetric because of the nonlocal condition (2), so, the structure of its spectrum is more complicated than for the Dirichlet conditions.The nonlocal integral conditions with γ never cause more problems than the classical conditions in both the number of the iterations and the precision of the solution.But these conditions affect the region of the convergence of the method.The convergence domain depends essentially on the coefficients of the nonlocality.The value of the parameter γ in nonlocal boundary conditions is essential for the stability of iterative methods.The results of the numerical experiment are in good agreement with the existing theoretical results for a two-dimensional elliptic equation in a rectangle domain with integral boundary conditions in one coordinate direction.

Lemma 11 .
If µ 1 > |η 1 |, then all eigenvalues of the matrix A 1 are positive.If µ 1 = |η 1 |, then one eigenvalue of the matrix A 1 is zero, other eigenvalues are positive.
2-order vector.We can verify directly that the next statements are valid.

Table 1 .
Accuracy of the solution and the number of the iterations for the case in which λ j (−Λx) > 0 for all j.Nonlinear Anal.Model.Control, 2014, Vol. 19, No. 3, 517-535

Table 2 .
The solution accuracy and the number of the iterations for the presence of at least one negative eigenvalue of the matrix −Λx: λ 1 (−Λx) < 0.