On the numerical solution for nonlinear elliptic equations with variable weight coefficients in an integral boundary conditions

Regimantas Čiupaila , Kristina Pupalaigė , Mifodijus Sapagovas Vilnius Gediminas Technical University, Saulėtekio ave. 11, LT-10223 Vilnius Lithuania regimantas.ciupaila@vilniustech.lt Department of Applied Mathematics, Kaunas University of Technology, Studentų str. 50, LT-51368 Kaunas, Lithuania kristina.pupalaige@ktu.lt Faculty of Mathematics and Informatics, Vilnius University, Akademijos str. 4, LT-08412 Vilnius, Lithuania mifodijus.sapagovas@mif.vu.lt

The boundary value problems for elliptic equations with nonlocal conditions as some elementary generalization of classical boundary value problems were formulated in [5,8].
On the other hand, many papers were published on the nonlocal boundary value problems for various types of equations as mathematical models of some physical phenomena related to plasma physics, heat transfer, thermoelasticity, chemical diffusion, underground water flow, biochemistry, population dynamics, etc. (see, e.g., [12] and references therein).
Both facts motivated active investigation of numerical methods for elliptic equations with nonlocal conditions. As a consequence, the numerical methods for linear and nonlinear elliptic equations with various nonlocal conditions were intensively investigated during past two decades. The numerical methods for linear two-dimensional elliptic equation with Bitsadze-Samarskii and multipoints nonlocal conditions were investigated in the papers [13,28,32].
The existence and uniqueness of a solution in the corresponding Sobolev spaces for a multidimensional elliptic equation with two integral boundary conditions ξ1 0 u(x, y) dy = 0, 1 ξ2 u(x, y) dy = 0, x = (x 1 , x 2 , . . . , x n ), were proved in [4]. The radial basis function collocation technique for the solution of such problem was presented in [18]. The convergence of finite difference method for the two-dimensional elliptic equation with condition (5) was proved in [7]. The nonlocal boundary conditions (2), (3) for the first time were formulated in [10,11] for the parabolic equation with application to thermoelasticity and thermodynamics.
The integral conditions in the form of (2), (3) are used by many authors considering nonlocal problems, also for differential equations of other types: for hyperbolic equations [16], equations with fractional derivatives [17], problems with complex parameter in the equation or nonlocal condition [23].
Essentially, condition (7) is a discrete analogue of condition (3). Considering difference schemes for the elliptic equation with nonlocal condition (7), many authors use the same assumptions for the coefficients α k as for the functions α(x), β(x) in integral conditions (2) and (3).
In the paper [1] the well-posedness of the second order of accuracy difference scheme for the elliptic equation is established under the assumption m k=1 |α k | < 1.
Analogous result for an inverse problem of multidimensional elliptic equation with nonlocal condition (7) is proved in [3] under the assumption m k=1 α k = 1, α k 0.
In [6] the slightly modified assumption which is sometimes more effective than other analogous assumptions, is used.
is presented in [14]. This assumption could be interpreted in a following way. If all the coefficients are nonnegative (α k 0), then (11) and (8) coincide. But if a part of the coefficients are negative (α k < 0), there is no limitation for them, and the sum of positive coefficients must not exceed one. With this assumption, the convergence of the difference scheme of the second order of accuracy is proved.
In the paper [29], several iterative methods for the solution of the system of difference equations approximating problem (1)-(4) are presented when Also we note that assumptions (6), (12) for the functions α(x) and ρ(λ) and assumptions (8)- (11) for the coefficients α k are only sufficient, but not the necessary conditions for the statements proved in the papers would be right.
The convergence of the difference schemes is a very important issue for any differential equation with nonlocal conditions. However, in the case of nonlinear elliptic equations, there is one more important task. It is necessary to consider how to solve the system of nonlinear difference equations. Usually, some iterative methods must be applied for this purpose. However, there are comparably not many papers in which the iterative methods for the difference systems with nonlocal conditions were investigated (see, e.g., [22,29] and references therein).
In [25] the Poisson equation with the variable α(x) and β(x) expressions in conditions (2), (3) is solved by the Peaceman-Rachford alternating direction method. The functions α(x) and β(x) are selected in a way that all the eigenvalues of respective difference operators would be positive. That guaranties the convergence of the method.
At the present time, it is not clear, where the assumptions, under which the convergence of difference schemes can be proved, are sufficient also for the convergence of some iterative methods. As far as authors know, that problem has not been investigated earlier.
It is important to emphasize that in most cases of elliptic equations with nonlocal conditions in the form of (2), (3) the matrix of difference problem under some assumption to functions α(x) and β(x) has the properties appropriate for the M-matrices [22,29]. Based on this property, it is possible to prove the convergence of many iterative methods [22].
Further, during some last years, the stability and convergence of difference schemes for some parabolic and elliptic equations with nonlocal conditions were investigated using the properties of the M-matrices [9,15,27]. We note that applications of the M-matrices for the elliptic and parabolic equations with Dirichlet boundary conditions have been described by Varga [30].
So, it is possible to state that the M-matrices theory is the appropriate methodology for other theoretical aspects of difference scheme in the case of nonlocal conditions, in particular, conditions in the form (2), (3).
The main aim of the present paper is to investigate when the matrix of the difference problem approximating problem (1)-(4) is an M-matrix with the weaker or more general assumptions for the coefficients α(x) and β(x) in comparison with the known assumptions (9)- (12). In the paper, using the properties of M-matrices, the new algorithm for the solution of the system of difference equations is provided.
The structure of the paper is following. In Section 2 the difference problem is formulated and rewritten in matrix form. The corresponding difference eigenvalue problem is considered in Section 3. Taking into apart, in this section the conditions, under which for all eigenvalues of difference operator the condition Re λ > 0 is satisfied, are obtained. In Section 4 the iterative methods for the solution of the system of nonlinear difference equations with nonlocal conditions are investigated. Results of numerical experiment supplementing the theoretical investigation are provided in Section 5. In Section 6 the generalizations and conclusions are formulated.

A difference problem
The boundary value problem with nonlocal conditions (1)-(4) is solved by the finite difference method. By this aim we provide that the unique, sufficiently smooth solution of this problem exists.
We rewrite the system of difference equations (19) in matrix form where A is matrix of the order (N − 1) 2 , f (u) and ϕ are vectors of the order (N − 1) 2 .
Vector f (u) is formed from the components f ij (u ij ), and vector ϕ-from the values of µ 1j ,μ 2j , µ 3i and µ 4i . We do not express the value of ϕ because it will not be used further.
Matrix A is formed as follows: where Λ is the square matrix of order (N − 1) 2 corresponding to a difference operator −δ 2 x − δ 2 y in the rectangular domain with the Dirichlet-type homogeneous boundary value conditions. Matrix C consists of multipliers hα i or hβ i of unknowns u ij in the equations of system (19). More exactly, C is a block matrix where The number of blocks of matrix C and order of matrix C 1 are N − 1.
Matrix A = Λ − C, taking into account the expressions of Λ and C, could be written also in other form: where I is an identity matrix of order (N − 1), https://www.journals.vu.lt/nonlinear-analysis Now we consider the eigenvalue problem for the matrix A that is necessary for the further investigation of the systems of nonlinear equations (19) or (20).

An eigenvalue problem of matrix A
We start from the difference eigenvalue problem We take this problem in matrix form. With this purpose, we express u 0j and u N j from conditions (26), (27): whereα i andβ i are defined by formulas (18). Putting these expressions of u 0j and u N j into equation (25) as i = 1 or i = N − 1, we get (24), where A is defined by formula (22). So, we get an important conclusion: Corollary 1. Difference eigenvalue problem (25)-(28) is equivalent to the eigenvalue problem for the matrix A.
Using the Fourier method, we separate variables in (25)-(28) Putting the expression of these unknowns into (25)- (28), analogously as in [24], we obtain two one-dimensional problems: and δ 2 y w j + ηw j = 0, j = 1, 2, . . . , N − 1, For the eigenvalues λ of problem (25)-(28), the following equality is true: For the corresponding eigenvectors, we have where i and j are the numbers of vectors coordinates, and k and l are numbers of eigenvectors.
The solution of eigenvalue problem (31) is known [19]: We rewrite problem (30) in matrix form analogously for problem (25)- (28). We express from nonlocal conditions of problem (30) whereα i andβ i are defined, as earlier, by formulas (18). Substituting these expressions into first equation of (30) as i = 1 and i = N − 1, we get where Λ x is defined by formula (23), Together with eigenvalue problem (30) or (33), we consider another problem as This problem could be written in matrix form as follows: We define auxiliary block-tridiagonal matrix Now we need a few properties of M-matrices. It follows from the definition that a kk > 0. We will denote A > 0 (A 0) if a kl > 0 (a kl 0). Additionally, A < B if a kl < b kl . We will use similar notation for vectors also.
If A is such that a kl 0, k = l, then the following statements are equivalent [31]: (S1) A −1 exists and A −1 0; (S2) The real parts of all the eigenvalues of the matrix A are positive: Re λ(A) > 0.
We will prove some properties of the matrices A and Λ x defined by formulas (22), (23).

Lemma 2. If hypotheses (H2) and (H3) are true and
then matrix Λ x defined by formula (23) is an M-matrix.
In this way the diagonal elements of matrix Λ (0) x are positive, and nondiagonal elementsnonpositive. It has been proved [20,21] that all the eigenvalues of the difference problem (34) are positive if and only if γ 1 + γ 2 < 2. So, according to Definition 1 and statements (S1), (S2), matrix Λ (0) x is an M-matrix. Further, we compare the elements of the matrices Λ x and Λ (0) x . From (18), according to hypothesis (H2), we get . Now it follows from (36) thatα α i . Analogously,β β i .
We denote elements of matrix Λ (0) x as b ij and elements of matrix Λ x as b ij . So, we have: x is an M-matrix.
It follows from these three properties that Λ x also is an M-matrix.
Remark 2. Condition (37) is only sufficient, but not necessary one for the matrix Λ x to be an M-matrix. This condition could be weakened with the concrete expressions α(x) and β(x). But in the general case, not concreting the expressions α(x) and β(x), is not possible to improve this condition. Indeed, e.g., when α(x) = M 1 , β(x) = M 2 , then Λ x is an M-matrix if and only if M 1 0, M 2 0, M 1 + M 2 < 2. And this is very close to condition (37).
We will obtain the conditions, which should be satisfied by the functions α(x) and β(x) that the eigenvalues of two-dimensional problem (25)-(28) satisfy the inequality Re λ kl > 0.
Proof of lemma is the same as for Lemma 3 in the paper [22] in which the case of γ 2 = 0, 0 γ 1 < γ 0 was investigated. There it was proved that where γ 0 is defined by formula (39), are true, then matrix A defined by (22) is an M-matrix.
Proof. Like in the proof of Lemma 2, γ 1 and γ 2 are defined by formulas (38). Then according to Lemma 3, we get λ kl (A (0) ) > 0 as condition (40) is satisfied. So, matrix Now we compare the elements of matrices A = {a ij } and A (0) . As for the elements of matrices Λ x = {b ij } and Λ (0) Remark 3. Assumption (40) is only sufficient but not necessary condition for the matrix A to be an M-matrix. If then λ(A (0) ) = 0 in the case (38). In this case, we could not exploit the method used in the proof of Theorem 1 for the investigation that the matrix A is an M-matrix or not. In more detail, this case will be investigated in Section 5.1 using the numerical experiment.
It is possible to reformulate the statement of Theorem 1. Precisely, the following statement is right.

Iterative methods
We solve the system of nonlinear difference equations (20) Au + f (u) = ϕ, where A = Λ−C, by iterative methods. According to Theorem 1, matrix A of this system defined by formula (22) is an M-matrix. We remind that for the nonlinear function f (u), hypothesis (H1) is satisfied. Proof. Statement of lemma follows from [22]. Now we can present the main result of our paper. In paper [22], under conditions α(x) 0, 1 0 α(x) dx ρ < 1, β(x) = 0, several implicit iterative methods for system (20) were investigated: where m 1 is constant from (H1). The main assumption for the iterative method to converge is following:  With these assumptions, matrix A is diagonally predominant. Meanwhile, with assumption (40), neither the symmetricity nor diagonal dominance are not necessary.
Each of these iterative methods is characterized by specific interpretation. To obtain u n+1 when the values of u n are known, we need to solve an additional simpler system of equations. The inner iteration is usually used for this. In method (41) the values u n+1 are obtained by solving the system of nonlinear difference equation with Dirichlet boundary condition. In method (42) the system of linear equation with nonlocal condition is solved. And using method (43), on every step of iteration the system of linear equation with Dirichlet boundary condition is solved. To obtain u n+1 when u n is known, in all three cases of iterative methods, it is necessary to solve the system of difference equations corresponding to two-dimensional elliptic equation.
Further, we consider one more iterative method in each step of which it is needed to solve only the system of linear tridiagonal equations with nonlocal conditions. With this aim, we write matrix A defined by formula (22) in the form where M and N are block-diagonal and block-tridiagonal matrices, respectively: We note that M is an M-matrix, N 0. Now for the system of equations (20), we write down the following iterative method: The convergence of method (44) is proved likewise as of iterative methods (41)-(43) in [22], on the basis of the properties of M-matrices. Iterative method (44) like methods (41), (42) and (43) is implicit one. The main advantage of method (44) is that on every step of the iteration, we have to solve not a twodimensional, but one-dimensional problem. We write method (44) by the coordinate form such as system (13)-(16):

Numerical results
In this section, we present some results of numerical experiment. The main goal of the numerical experiment is to illustrate the theoretical results obtained in the paper. Our numerical experiment consists of two parts.

Conditions under which the matrix of difference problem is an M-matrix
This is the first part of numerical experiment. In Lemma 2 and Theorem 1 the sufficient conditions (37) and (40) for the matrix Λ x and respectively A were formulated under which these matrices are the M-matrices. Condition (40) differs from the ones used by other authors investigating difference schemes for two-dimensional differential equations with nonlocal boundary conditions of the same form (see (6), (8)- (12)).
To explain the effectiveness of assumption (40), we carried out the numerical experiment with fixed values of h and some different expressions of α(x) and β(x) directly computing all the eigenvalues of matrix Λ x .
To obtain the eigenvalues of matrix A, we used formula (32). For chosen expressions α(x) and β(x), we varied coefficients M 1 and M 2 in these expressions checking when property is satisfied. In this way, using the numerical experiment, we obtain the area on the coordinate plane (M 1 , M 2 ), where the inequality Re λ(A) > 0 is satisfied for all the eigenvalues of A. We compared the results obtained with the theoretical results presented in Theorem 1. We chose five following different functions α(x) and β(x).  according to the results of the numerical experiment. It means that the matrix A in this (bigger) area is an M-matrix according to the results of the numerical experiment.
The following conclusions could be formulated from results of the numerical experiment.
If the functions α(x) and β(x) change sufficiently slow in the interval x ∈ (0, 1), i.e., when the variations of these functions are relatively insignificant, then the sufficient condition (40) for A to be an M-matrix is close to be the necessary condition. In other words, the pointwise curve in this case is close to the continuous curve. This conclusion is illustrated by the results of numerical experiment provided in Figs. 1(a) and 1(d) for Examples 1 and 4. In Example 1 the variation of functions α(x) and β(x) is half of size of these functions in Examples 2, 3, 5.
If functions α(x) and β(x) sufficiently differ from constants, assumption (40) remains correct, but it could be noticeably improved. From Figs. 1(b), 1(c), 1(e) it could be seen that matrix A is an M-matrix (Re λ(A) > 0) with less restricted assumption for M 1 and M 2 in comparison with the statement of Theorem 1.

Iterative methods
In the second part of the numerical experiment the system of nonlinear difference equations (13)-(16) was solved by iterative method (45).
The values of M 1 and M 2 in the expressions of functions α(x) and β(x) were chosen from the first part of numerical experiment. More precisely, the theory (Theorem 1) does not guarantee that with the chosen functions α(x) and β(x), matrix A is an M-matrix. However, according to the results of numerical experiment, it distinguishes by such property. Thus, iterative method (44) must converge, what is illustrated by the results of computations.
We consider problem (1)-(4) with f (x, y, u) = 2u 3 e −2(x+y−1) , These expressions were chosen in a way that the function u(x, y) = e x+y−1 would be the exact solution of this problem. We admit that functions α(x) and β(x) are chosen in the same way as in Example 5.
Numerical results are presented in Tables 1-3.  In all Tables 1-3 the errors = max (i,j) |u ij −u * (x i , y j )| of the approximate solution u ij are presented, where u * (x i , y j ) is the solution of the differential problem.
The results provided in Table 1 show the dependence of the error from α(x) and β(x) when M 1 and M 2 satisfy the condition M 1 + M 2 = 3.42.   The iterative method converges, the error is of the order O(h 2 ). We admit that the error is decreasing lawfully when M 1 is increasing (correspondingly, M 2 is decreasing). In other words, inequality M 1 > M 2 influences smaller error in comparison with the case of M 1 < M 2 .
In Table 2 the values of the errors are provided for some qualitative different values of (M 1 , M 2 ). Namely, if M 1 = M 2 = 0.1 then nonlocal conditions (2), (3) are close to Dirichlet conditions. If M 1 = 2, M 2 = 1.42, the matrix A is an M-matrix according to Theorem 1. And in the case of M 1 = 0.42, M 2 = 10.1, condition (40) of Theorem 1 is not satisfied, but the matrix A is an M-matrix according to the results of the numerical experiment. Again, we admit that the error is of the order O(h 2 ), i.e., the error is proportional to the error of approximation.
In Table 3 the results on the dependence of the error on the values M 1 and M 2 are provided when M 1 = M 2 . This error is steadily, although slowly increasing when the value M 1 = M 2 is increasing. This could be explained by the fact that the error of trapezoid formula increases when M 1 = M 2 increases.
We could also formulate one more subtle conclusion. When the point (M 1 , M 2 ), M 1 = M 2 is placed above of the pointwise line, matrix A is no longer an M-matrix according to the results of numerical experiment. Therefore, condition (40) is not satisfied. However, the iterative method still converges for some value M 1 = M 2 a little bit bigger than 3.42. It could be explained because the convergence of the iterative method (44) is assured by the matrix M + m 1 I is being an M-matrix.