M-matrices and Convergence of Finite Difference Scheme for Parabolic Equation with an Integral Boundary Condition

In the paper, the stability and convergence of difference schemes approximating semilinear parabolic equation with a nonlocal condition are considered. The proof is based on the properties of M-matrices, not requiring the symmetry or diagonal predominance of difference problem. The main presumption is that all the eigenvalues of the corresponding difference problem with nonlocal conditions are positive.


Introduction
Initiations in the research of parabolic equations with nonlocal conditions appeared in papers [3,19] more than half a century ago. New problems, formulated in the papers mentioned, became as the object of investigations of many authors later on.
During the last decades many new papers have appeared, where solution of parabolic equations with integral boundary conditions by the finite difference method was considered in various aspects, including the stability and convergence of difference schemes. In paper [11], the estimation of error in the maximum norm was obtained for the one-dimensional linear parabolic equation with integral boundary conditions.The analogous result was provided in [12] for the numerical solution of a semilinear problem with integral boundary conditions. Various aspects of difference schemes, such as convergence, stability, algorithms of realization, approximation of increased accuracy in one-dimensional case were investigated in many articles (see, for example [8,9,22,26] and the references therein). The influence of complex coefficients on the stability of difference schemes was investigated in [20,31].
In a two-dimensional case of parabolic equations with integral conditions, the difference methods were applied in [4,6,21,24,29], and to the equations with other type of nonlocal conditions -in [13,15].
We emphasize the fact that most of theoretical issues of difference schemes were investigated in line with the sufficiently strong conditions introduced (it is not clear whether they are necessary) for the various parameters or functions in nonlocal conditions (see, for example [11,12,21]), or proving the stability or convergence in unusual energetic norms (see [13,16]).
Another important aspect is that theoretical investigations in the numerical analysis are strongly influenced by practical applications in heat conductions [3], thermoelasticity [10], underground water flow [23], electrochemistry [5], and so on. In [7,14], some mathematical models with nonlocal conditions for bioreactors are described.
In this paper, the stability and convergence of difference schemes for twodimensional nonlinear parabolic equations with an integral boundary condition are proved using the properties of M-matrices and the structure of spectrum of difference operators with nonlocal conditions.
In papers [28,33,38], the theory of M-matrices was used for the investigation of difference systems obtained from the elliptic equation with an integral condition. In many cases, the matrix of the system of difference equations, approximating a differential equation with nonlocal conditions, is characterized by the properties specific to M-matrices. These properties might be used to obtain the conditions of convergence of iterative methods [28,38]. In the paper [18], the properties of M-matrices were applied to the investigation of the stability of difference schemes in the energetic norm.
In the paper [33], using the properties of M-matrices, the error of the solution of the difference problem for nonlinear elliptic equation with nonlocal boundary condition was estimated. Two aspects of such a method of proof have been noted. First, the main idea of error estimation is the construction of majorant, similar as applying the maximum principle method. Second, in the M-matrix method, the error estimate is proved declining the diagonal predominance property of the matrix for the system of difference equations, which is necessary for the maximum principle method. The results of the investigation of the structure of the difference operator spectrum are used instead of the diagonal predominance property.
The stucture of the spectrum of difference operators with nonlocal conditions was investigated in many papers (see [30,35,36] and also references in review article [37]). In the papers [2,16,17,27,34], the stability of difference schemes in the energetic norm according to the structure of the spectrum is examined.
In this article the methodology of the estimation of the error of the solution as in [33] is applied to a new class of problems with nonlocal conditions, namely, to nonlinear two-dimensional parabolic equations. The main result of the present paper is that the stability and convergence of difference schemes for a parabolic equation with an integral boundary condition in the uniform norm has been proved using the structure of the spectrum of difference operator and the theory of M-matrices. The statement on the convergence of the difference scheme in the uniform norm was proven only in the case of nonlocal conditions under which the matrix of the system of difference equations is diagonally dominant (see, for example [11,12,21]). Using the methodology of M-matices the requirement of diagonal dominance was abolished.
The paper is organized as follows. In Section 2, the differential and difference problems are stated. The difference problem for the error of the solution is investigated in Section 3. The latter problem is expressed in the matrix form as a two-layer difference scheme. In Section 4, some new formulations of the well-known properties of M-matrices, adapted to the system of difference equations are obtained. Using these properties, the error estimation is presented for the solution and the convergence of the difference scheme in the uniform norm is proved in Section 5 . In Section 6, a short analysis of the stability of a difference scheme is provided. Results of numerical experiments are delivered in Section 7. Some final conclusions are presented in Section 8.
The following notation is presented A finite difference scheme is now defined Further in this article for shorter writing no restrictions are imposed on index n in difference equations or conditions, assuming that n is any integer 1 ≤ n ≤ N 1 . The solution U n−1 i,j at the all points (i, j) of the layer t = t n−1 is known. Then system (2.5)-(2.8) for the unknowns U n ij with fixed n can be interpreted formally as difference analogue of some elliptic equation. The solution of such system by the iterative methods was investigated in [28]. It follows from these results that under hypotheses H1 and H2, the unique solution U n ij of the system of difference equations (2.5)-(2.8) exists.
Thus the following statement is obtained.
Lemma 1. If the hypotheses H1 and H2 are true, then unique solution of the system of difference equations (2.5)-(2.8) exists.
Let the differential problem (2.1)-(2.4) possess the unique smooth enough solution, i.e. the derivatives of the solution with respect to x and y up to the fourth order and the first and second derivatives with respect to t are bounded. Then the error of approximation is O(h 2 + τ ). So, the following difference problem for the solution u n ij of the differential problem can be written: Here, according to the presumption on the smoothness of solution of the differential problem, the following estimations are true:

Difference problem for the error
The error of the solution is noted z n ij = u n ij −U n ij . From the problems (2.5)-(2.8) and (2.9)-(2.12) the following system recieved: where the following is denoted andũ n ij is a certain intermediate valueũ n ij = u n ij + θU n ij , |θ| ≤ 1. To estimate the error z n ij , first of all this system is rearranged. In this system, where n is a fixed number, there are N (N − 1) equations and the same number of unknowns. This system is reduced to two separate systems of lower order. With this aim z n 0j is expressed from equation (3.2) for each j = 1, 2, . . . , N − 1: As h ≤ 1 2 and according to hypothesis Putting expression (3.6) into equations (3.1), as i = 1 and introducing the new notation, for each fixed n ≥ 1 is obtained:  -(N − 1) 2 order system (3.8), (3.9), (3.4) with the unknowns z n ij in the interior points of the grid area, i.e. as i, j = 1, 2, . . . , N − 1 and -N − 1 order system (3.6) with the unknowns z n 0j , i, j = 1, 2, . . . , N − 1; indeed, this system is the explicit formulas providing the opportunity to calculate z n 0j , when the solution of the first system z n ij , i, j = 1, 2, . . . , N − 1 is known.
Now the system (3.1), (3.2), (3.4) for fixed n ≥ 1 is written in the matrix form Az n = Bz n−1 + r n , Matrices I, Λ, C and D of order (N − 1) 2 included to this system are formed as follow. I is the unique matrix. Λ is the matrix, corresponding to the difference Laplace operator (with a minus sign) with zero Dirichlet conditions. C is the matrix composed of the coefficients in the member h −2 α z n ij of equation (3.8). D is the diagonal matrix with the elements d n ij , prescribed by formula (3.5). For details see [33]. Denote Thus, the system of difference equations is written in the form of (3.10). Another form of this system is  In the present paper, the following notation is used: A > 0 (A ≥ 0) if a kl > 0 (a kl ≥ 0) for all k, l. Also A < B (A ≤ B), if a kl < b kl (a kl ≤ b kl ) is written. Analogous notation will also be used for the vectors.
The following three lemmas could be interpreted as the new statement of certain properties of M-matrices, adapted to the systems of difference equations approximating parabolic equations. As far as the authors are acquainted, these properties of M-matrices were not formulated in this form so far.

Lemma 2.
Let v n be the solution of the difference equation Av n = Bv n−1 + f n , n ≥ 1. (4.1) If A is an M-matrix, B ≥ 0, v 0 ≥ 0 and f n ≥ 0 for all n ≥ 1, then v n ≥ 0 for all n ≥ 1.
Proof. The statement of Lemma 2 follows by using the method of mathematical induction. We Suppose there is the vector v with the coordinates v k , k = 1, 2, . . . , m. Vector with the coordinates |v k | is designated by |v|, i.e. |v| = {|v k |}. respectively. If A is an M-matrix, B ≥ 0, w 0 ≥ 0, g n ≥ 0 as n ≥ 1 and, additionally, |v 0 | ≤ w 0 , |f n | ≤ g n , n ≥ 1, then |v n | ≤ w n , n ≥ 1.
Proof. When D ≥ 0 is a diagonal matrix and A is an M-matrix, then A + D is also an M-matrix. Moreover Similary as in the proofs of Lemmas 2 and 3, from |v n−1 | ≤ w n−1 it follows

Error estimation and convergence of the difference scheme
With reference to the results of Section 4 , we estimate the solution of the system of difference equations (3.1)-(3.4), i.e. the error z n ij = u n ij − U n ij . In Section 3, for each fixed n ≥ 1 this system was reduced to two systems: system (3.2), (3.9), (3.4) and explicit formulas (3.6). The first system was presented in matrix form (3.10), where matrices A and B were described by formulas (3.11), and the vectorr n was composed of the expressions of the right-hand sides of equations (3.4), (3.9), i.e.r n = {r n ij }: Proof. In the article [33] is proved, that the matrix Λ − C under conditions H1 and H2 are M-matrices (see [33] Corollary 1). Since I ≥ 0 and D ≥ 0 are diagonal matrices, it follows from the properties of M-matrices that matrices A 1 = τ −1 I + Λ − C and A = τ −1 I + Λ − C + D are M-matrices as well. Proof. First of all, it is estimated that the solution z n ij of this system, as i, j = 1, 2, . . . , N − 1. In other words, in the beginning, the solution of system (3.8), (3.9), (3.4) is estimated. Two separate cases: 0 ≤ γ ≤ 1 and 1 ≤ γ ≤ 2−ρ are considered.
Thus, for the systems of equations (3.10) and (5.14) all the presumptions of Lemmas 3 and 4 are correct. Hence it is derived |z n ij | ≤ w n ij , i, j = 1, 2, . . . , N − 1, n = 1, . . . , N 1 . In such way, from the expression of the function w(x, y, t) the following recieved |z n ij | ≤ Now, according to formula (3.6), it is found 6 The stability of a difference scheme Estimation (5.17) of error could also be interpreted as statement on stability of the difference scheme. Explaining more in detail, the definition of stability is introduced. Since difference scheme (2.5)-(2.8) is nonlinear, the common concept of stability can be used. Namely, two problems are considered. The first problem is (2.5)-(2.8) with the data p n ij , (µ k ) n j , k = 1, 2, (µ l ) n i , l = 3, 4 and ϕ ij . The solution of this problem is U n ij . The second problem is taken with the modified datap n ij , (μ k ) n j , (μ l ) n i andφ ij . The solution of this problem is denoted byŨ n ij . Definition 2. Difference scheme (2.5)-(2.8) is stable if for every real ε > 0 the value δ = δ(ε), not depending on h and τ , exists such that |U n ij −Ũ n ij | ≤ ε, i, j = 0, 1, . . . , N ; n = 1, 2, . . . , N 1 if the data of both problems considered differs by not more than δ. Proof. According to Definition 2, it is necessary to prove that for each ε > 0 the value δ = δ(ε), not depending on h and τ , will arise, such that, for the solution of system of difference equations (6.1)-(6.4), the estimate is correct if (6.5) is valid. This statement could be fully proved analogously as Theorem 1. It is only enough to change the definition of the function w(x, y, t) a little. Namely in the definition of the function w 1 (x, t) (formula (5.2)) the constant δ can be used instead of C2τ 2 . In the definition of the function w 2 (x, y) (formula (5.3)) can be used δ instead of constant M h 2 24 . Then analogously to estimates (5.16) and (5.17), it is obtained |z n ij | ≤ (T + 1)δ + 13δ/ρ = (T + 1 + 13/ρ) δ, i, j = 1, 2, . . . , N − 1 (6.7) and |z n 0j | ≤ (8T + 12 + 13/ρ) δ, j = 1, 2, . . . , N − 1.

Numerical results
In order to verify the theoretical results and demonstrate the order of convergence of the approximate solution, some test problem where the exact solution of differential problem is known is considered. A differential equation is considered 3), boundary conditions (2.4) and initial condition (2.2). The functions p(x, y, t), ϕ(x, y), µ i (y, t), i = 1, 2, and ν i (x, t), i = 3.4 are chosen that the function u(x, y, t) = e at sin πx sin πy + bx 2 (7.1) is the analytical solution of the problem under investigation. The test problem has been solved with some different parameters γ, a, b, T , h and τ . In order to evaluate the accuracy of the numerical method, the absolute or relative error were used The theoretical results presented in the previous section do not depend upon the numerical method that is used to solve system of nonlinear difference equations (2.5)-(2.8) with fixed n.
This system can be solved using one of the iterative methods suitable for problems with nonlocal conditions [28]. Here a Jacobi method was used.
The results of the numerical test for different h and τ are presented in Table 1. It reveals that the absolute error E n h,τ decreases approximately as O(h 2 + τ ). Now the role of additional term bx 2 in (7.1) is explained. The  Table 2. In Table 3, the results of numerical experiment of influence of T on the error of solution is presented. All constants M 2 , M 4 and C 2 from (2.13) depend on e at for solution (7.1). It means that in the case a > 0 these constants as well as the absolute error E n h,τ for fixed h and τ can grow rapidly, when t is growing. In case a > 0 the relative error ε n h,τ is used. The relative error with fixed h and τ increase slowly with the growing t (see the estimate of constant C 4 in (5.17), which depends on T ). The numerical results confirm the stability of the method.

Conclusions
The results obtained in the paper extend the class of the differential problems with nonlocal conditions when it is possible to prove the stability and convergence of difference schemes in the uniform norm. In this case, there is no need to require the matrix of the difference problem to be symmetrical or diagonally dominant. In this paper, the conditions of the stability and convergence are received using the main presumption that all the eigenvalues of the corresponding difference problem with nonlocal conditions are positive. Thus, nonlocal condition (2.3) can be interpreted as one particular case of nonlocal conditions when it is possible to prove the stability and convergence of the difference scheme in the uniform norm. So, the methodology of M-matrices might be also used in different cases of nonlocal conditions when all eigenvalues of the difference problem are positive. Some of such cases are examined in review article [37].
Furthermore, from the properties og M-matrices it follows that it would be enough to require even more general conditions. Namely, that the condition Reλ > 0 should be true for all the eigenvalues of the corresponding matrix. In papers [30,32,34], several specified examples with other types of nonlocal conditions are provided, when there exist complex eigenvalues, the real parts of which are only positive. It might be possible to apply the methodology of the present paper to these problems. The relevant majorant w(x, y, t) should be constructed.
Regarding, that estimates (5.17), and (6.7), (6.8), obtained in the present paper, depend on the constant ρ (see the hypothesis H2 in Section 2). It is not clear whether this constant should be certainly included in the estimate. It is only possible to state that the estimates depend on the choice of the majorant w(x, y, t). In addition, the function w 2 (x, y) is very close to the majorant used in the maximum principle to the stationary problems [25].