Electronic Journal of Differential Equations, Vol. 2022 (2022), No. 44, pp. 1–15. ISSN: 1072-6691. URL: http://ejde.math.txstate.edu or http://ejde.math.unt.edu STABILITY ANALYSIS OF THE PEACEMAN-RACHFORD METHOD FOR PARABOLIC EQUATIONS WITH NONLOCAL CONDITIONS MIFODIJUS SAPAGOVAS, JURIJ NOVICKIJ, REGIMANTAS ČIUPAILA Abstract. We consider an efficient finite difference method solving of two- dimensional parabolic equations with nonlocal conditions. The specific feature of the investigated problem is that the nonlocal condition contains the values of solution’s derivatives at different points. We prove the stability of this method in specific energy norm. The main stability condition is that all eigenvalues of the corresponding difference problem are positive. Results of computational experiments are presented. 1. Introduction In this article, we solve a two-dimensional linear parabolic equation using finite difference method (FDM), ∂u ∂t = ∂2u ∂x2 + ∂2u ∂y2 + f(x, y, t), (x, y) ∈ Ω, t ∈ (0, T ], (1.1) where Ω = (0, 1)× (0, 1), with nonlocal boundary conditions γ ∂u(0, y, t) ∂x = ∂u(1, y, t) ∂x , γ ∈ (0, 1], (1.2) u(x, 0, t) = u(x, 1, t) = u(0, y, t) = 0, (1.3) and initial condition u(x, y, t) = ϕ(x, y), (x, y) ∈ Ω. (1.4) Intensive research of differential problems with various nonlocal conditions started after Cannon and Kamynin had published their works [5, 16]. The authors, instead of classical boundary condition u(0, t) = µ(t), formulate nonlocal condition∫ x 0 u(x, t)dx = E(t), where E(t) is the known amount of heat in interval (0, x) at time t. To the authors’ knowledge, boundary condition of type (1.2) first time was for- mulated in [14] for one-dimensional parabolic equation. It was observed that mathe- matical models with nonlocal conditions (1.2) are encountered in physical problems describing processes of particles’ diffusion in turbulent plasma. 2020 Mathematics Subject Classification. 65M06, 35K20, 34B10, 34K20. Key words and phrases. Nonlocal boundary conditions; parabolic equations; alternating direction method; stability of finite difference scheme. ©2022. This work is licensed under a CC BY 4.0 license. Submitted October 9, 2021. Published June 30, 2022. 1 2 M. SAPAGOVAS, J. NOVICKIJ, R. ČIUPAILA EJDE-2022/44 For the one-dimensional problem ∂u ∂t = ∂2u ∂x2 , 0 < x < 1, t > 0, u(x, 0) = u0(x), u(0, t) = 0, γ ∂u(0, t) ∂x = ∂u(1, t) ∂x , γ ∈ (0, 1], the stability of finite difference schemes (FDS) is proved in rather complex energy norms (see [9, 10, 12, 13] and references therein). The structure of FDS spectrum is investigated in detail for this purpose. The equivalence between these energy norms and L2-norm in vector space is proved in [10]. Subject of various type differential problems with nonlocal conditions like the uniqueness and existence of a solution, numerical methods and applications, spec- trum analysis is widely investigated (see e.g. [2, 8, 18, 19, 21, 33] and references therein). Problems with various types of nonlocal conditions are considered to be one of modern areas of research in theory of differential equations and numerical analy- sis. Mentioned research area is bound together with applications in science and technology. Nonlocal problems have special feature. The structure of spectrum of differential or difference operators with nonlocal conditions is more complicated and substantive compared to respective spectrum of operator with classical (Dirichlet or Neumann) boundary conditions (see e.g. [3, 17, 23, 28, 33] and references therein). Problem (1.1)–(1.4) with f = 0 is solved in [11] using Crank-Nicolson difference method with approximation error O(h2 + τ2). Stability of used method is proved in special energy norm ‖u‖D = (Du, u)1/2, where D is symmetric positive definite matrix depending on eigenvectors and asso- ciated vectors of difference problem. It is well known that stability is one of the most important properties of numerical method both for theoretical research and practical applications. Another important point arises for the two- and multi-dimensional parabolic problems with nonlocal conditions: how to solve finite difference scheme efficiently at every time layer tk, when the method is implicit. The main task of this article is to construct efficient FDM for the two-dimensional parabolic equation (1.1) with nonlocal boundary condition (1.2). The example of FDM for problems with classical (Dirichlet or Neumann) boundary conditions is Peaceman-Rachford alternating direcion implicit (ADI) method. In this article we modify Peaceman-Rachford ADI method by applying it to the special nonlocal condition (1.2) and investigate stability conditions. For this pur- pose we use results about spectrum structure of one-dimensional difference problem with nonlocal conditions. To the authors’ knowledge, ADI method for the parabolic equation with (1.2) type nonlocal condition has not been investigated earlier. For the two-dimensional parabolic and elliptic equations with other type nonlocal conditions ADI methods theoretically and practically were investigated in [7, 27, 30, 34]. This article is organized as follows. In Section 2 we formulate Peaceman-Rachford ADI method for differential problem (1.1)-(1.4). We prove that approximation er- ror is O(h2 + τ2). In Section 3 we investigate stability of the method in a special EJDE-2022/44 PEACEMAN-RACHFORD METHOD FOR NONLOCAL PROBLEMS 3 energy norm. We provide numerical examples for the ADI method in Section 4. Remarks and generalizations are provided in Section 5. 2. Numerical method and approximation error First, we define Unij := U(xi, yj , t n), where xi = ih, i = 0, N ; yj = jh, j = 0, N ; h = 1/N ; tn = nτ , n = 0,M ; τ = T/M ; N,M ∈ Z. We denote differences approximating derivatives of solution as δxU n ij := Uni+1,j − Unij h , δx̄U n ij := Unij − Uni−1,j h , δtU n ij := Un+1 ij − Unij h , δ2 xU n ij := Uni−1,j − 2Unij + Uni+1,j h2 , δ2 yU n ij := Uni,j−1 − 2Unij + Uni,j+1 h2 . We formulate Crank-Nicolson difference method for the problem (1.1)–(1.4). We emphasize that approximation accuracy of nonlocal condition (1.2) is O(h2 + τ2). Un+1 ij − Unij τ = 1 2 ( δ2 xU n+1 ij + δ2 yU n+1 ij ) + 1 2 ( δ2 xU n ij + δ2 yU n ij ) + f n+1/2 ij , i, j = 1, N − 1, (2.1) Un+1 Nj − UnNj τ = 1 2 ( 2 h ( γδxU n+1 0j − δx̄Un+1 Nj ) + δ2 yU n+1 Nj ) + 1 2 ( 2 h ( γδxU n 0j − δx̄UnNj ) + δ2 yU n Nj ) + γf n+1/2 0j + f n+1/2 Nj , j = 1, N − 1, (2.2) Un+1 0j = Un+1 i0 = Un+1 iN = 0, i, j = 0, N, (2.3) U0 ij = ϕij , i, j = 0, N. (2.4) Equation (2.1) is written in standard form changing second derivatives with differences. It approximates differential equation (1.1) at all inner Ω points i, j = 1, N − 1 with accuracy O(h2 + τ2), under the assumption that the solution of differential problem is sufficiently smooth. Equation (2.2), approximating nonlocal condition (1.2), is derived in the follow- ing way. First, condition (1.2) is rewritten as γ ( δxu n+1/2 0j − h 2 ∂2u n+1/2 0j ∂x2 +O(h2) ) = δx̄U n+1/2 Nj + h 2 ∂2u n+1/2 Nj ∂x2 +O(h2). (2.5) We make an assumption that differential equation (1.1) is valid not only in inner Ω points, but also on boundaries, when x = 0 and x = 1 (i.e. i = 0 and i = N). This assumption is usual for FDM if boundary condition has solution’s derivative with respect to spatial variable [24]. In accordance with this assumption, we re- place second order derivatives with respect to x in (2.5) by their expression from differential equation (1.1). Further step is approximation of differential expression: derivatives ∂u/∂t and ∂2u/∂y2 are replaced by differences with accuracy O(h2+τ2). We also use the fact that Uij n+1/2 = 1 2 ( Un+1 ij + Unij ) +O(τ2). We obtain equation 4 M. SAPAGOVAS, J. NOVICKIJ, R. ČIUPAILA EJDE-2022/44 an approximating nonlocal condition (1.2) with accuracy O(h2 + τ2) by eliminat- ing approximation errors. After some trivial transformations we obtain difference equation of form (2.2). Now we rewrite equations (2.1) and (2.2) in the form Un+1 ij − Unij τ = 1 2 ( δ̃2 xU n+1 ij + δ2 yU n+1 ij ) + 1 2 ( δ̃2 xU n ij + δ2 yU n ij ) + f̃ n+1/2 ij , (2.6) where δ̃2 xU n ij := { δ2 xU n ij , i = 1, N − 1, 2 h ( γδxU n 0j − δx̄UnNj ) , i = N, (2.7) f̃ n+1/2 ij := { f n+1/2 ij , i = 1, N − 1, γf n+1/2 0j + f n+1/2 Nj , i = N. (2.8) Remark 2.1. As noted above, difference equation (2.2) approximates nonlocal condition (1.2) with accuracy O(h2 + τ2), assuming equation (1.1) is valid and on boundary x = 0 and x = 1. The advantage of this approach is that difference equation (2.2) formally (but only formally) has the same structure, as difference equation (2.1). Actually, term δ̃2 xU n Nj is not an approximation of ∂2U/∂x2 at the point (xN , yj , t n). But, this term is a kind of generalized difference analogue of derivative. Expanding Un1j and UnN−1,j in the Taylor series, when |∂3u/∂x3| ≤ C and nonlocal condition (1.2) is valid, we have δ̃2 xU n Nj = 2 h ( γ Un1j − Un0j h − UnNj − UnN−1,j h ) = 2 h ( γ ∂Un0j ∂x + γ h 2 ∂2Un0j ∂x2 − ∂UnNj ∂x + h 2 ∂2UnNj ∂x2 +O(h2) ) = γ ∂2Un0j ∂x2 + ∂2UnNj ∂x2 +O(h). (2.9) Similarly, we can write another form of this expression δ̃2 xU n Nj = γ ∂2Ũn0j ∂x2 + ∂2ŪnNj ∂x2 , (2.10) where Ũn0j = U(x̃0, yj , t n), x̃0 ∈ [0, h]; ŪnNj = U(x̄N , yj , t n), x̄N ∈ [1 − h, 1]. This remark will be used for evaluation of ADI method’s aproximation error. Now, taking into account approximation form of nonlocal condition (1.2) in Crank-Nicolson method, we construct a following alternating direction method for differential problem (1.1)–(1.4) (Peaceman-Rachford ADI method) U n+1/2 ij − Unij τ/2 = δ2 xU n+1/2 ij + δ2 yU n ij + f n+1/2 ij , i, j = 1, N − 1, (2.11a) U n+1/2 Nj − UnNj τ/2 = 2 h ( γδxU n+1/2 0j − δx̄Un+1/2 Nj ) + δ2 yU n Nj + γf n+1/2 0j + f n+1/2 Nj , j = 1, N − 1, (2.11b) U n+1/2 0j = 0, j = 1, N − 1, (2.11c) Un0j = Uni0 = UniN = 0, (2.11d) EJDE-2022/44 PEACEMAN-RACHFORD METHOD FOR NONLOCAL PROBLEMS 5 Un+1 ij − Un+1/2 ij τ/2 = δ2 xU n+1/2 ij + δ2 yU n+1 ij + f n+1/2 ij , i, j = 1, N − 1, (2.12a) Un+1 Nj − U n+1/2 Nj τ/2 = 2 h ( γδxU n+1/2 0j − δx̄Un+1/2 Nj ) + δ2 yU n+1 Nj + γf n+1/2 0j + f n+1/2 Nj , j = 1, N − 1, (2.12b) Un+1 i0 = Un+1 iN = 0, j = 1, N − 1, (2.12c) U n+1/2 0j = U n+1/2 i0 = U n+1/2 iN = 0. (2.12d) The method, described for formulas (2.11) and (2.12), can be rewritten in another form taking into account (2.7) and (2.8). We will use the following form in further theoretical investigation U n+1/2 ij − Unij τ/2 = δ̃2 xU n+1/2 ij + δ2 yU n ij + f̃ n+1/2 ij , (2.13) Un+1 ij − Un+1/2 ij τ/2 = δ̃2 xU n+1/2 ij + δ2 yU n+1 ij + f̃ n+1/2 ij , (2.14) for i = 1, N , j = 1, N − 1. One needs to solve N -th order system N − 1 times (for every fixed j = 1, N − 1) with unknowns U n+1/2 1j , U n+1/2 2j , . . . , U n+1/2 Nj to find U n+1/2 ij from system (2.13). Similarly, one can find Un+1 ij from system (2.14). The matrix of system (2.14) is tridiagonal and boundary conditions are of Dirich- let type. The matrix of system (2.13) differs from diagonal matrix only by one element (U n+1/2 1j coefficient in the N − th equation is not equal to zero). Systems of such type usually arise in solving one-dimensional boundary value problem with periodic boundary condition [26]. The number of arithmetic operations to solve systems (2.13) and (2.14) is proportional to N . Therefore, number of arithmetic operations to find Un+1 ij , when Unij are known, using algorithm (2.13), (2.14) is pro- portional to N2 (that is proportional to the number of grid points). Such algorithms are called efficient (or economical, see e.g. [24]). Each one of the systems (2.13) and (2.14) separately approximates differential problem (1.1)–(1.4) with approximation error O(h2 + τ). Now, we prove that the system (2.13)–(2.14) approximates differential problem with accuracy O(h2 + τ2). For this purpose we use the method described in [24, Ch. 9, §1] for parabolic equations with boundary conditions of Dirichlet type. We eliminate intermediate values U n+1/2 ij from the system (2.13)–(2.14). For this purpose, we substract eq. (2.13) from (2.14) for every index pair (i, j). For every i value we have U n+1/2 ij = Un+1 ij + Unij 2 − τ 4 δ2 y ( Un+1 ij − Unij ) , i = 1, N, j = 1, N − 1. (2.15) 6 M. SAPAGOVAS, J. NOVICKIJ, R. ČIUPAILA EJDE-2022/44 By substituting this expression into system (2.13), after some transformations, we obtain Un+1 ij − Unij τ = 1 2 ( δ̃2 xU n+1 ij + δ2 yU n+1 ij ) + 1 2 ( δ̃2 xU n ij + δ2 yU n ij ) − τ 4 δ̃2 xδ 2 y Un+1 ij − Unij τ + f̃ n+1/2 ij , i = 1, N, j = 1, N − 1. (2.16) Systems (2.13)–(2.14) and (2.16) are equivalent. Now, we can evaluate approximation error of the ADI method (2.13)–(2.14). Suppose one more assumption is valid (except for the standard differential problem’s smoothness conditions, which ensure approximation order O(h2 + τ2) of Crank- Nicolson method) | ∂5u ∂x2∂y2∂t | ≤ C. Then, the corresponding differences are also bounded, and ∣∣δ2 xδ 2 y Un+1 ij − Unij τ ∣∣ ≤ C regardless of h and τ values.The following inequality is also valid (according to (2.10)) ∣∣δ̃2 xδ 2 y Un+1 ij − Unij τ ∣∣ ≤ 2C. Therefore, the ADI method written in the form of (2.16) differs from Crank- Nicolson method only in O(τ2) term. The approximation error of both methods is of the O(h2 + τ2) order. 3. Stability of difference scheme In this section we investigate the stability of the ADI method (2.13)–(2.14) with boundary conditions (2.11c), (2.12c). First, we rewrite the ADI method in the matrix form. We define N ×N matrix Λx and (N − 1)× (N − 1) matrix Λy by Λx =  2 −1 0 · · · 0 0 −1 2 −1 . . . 0 0 0 −1 2 . . . 0 0 ... . . . . . . . . . . . . ... 0 0 0 −1 2 −1 −2 0 0 · · · −1 2  , Λy =  2 −1 0 · · · 0 0 −1 2 −1 . . . 0 0 0 −1 2 . . . 0 0 ... . . . . . . . . . . . . ... 0 0 0 −1 2 −1 0 0 0 · · · −1 2  EJDE-2022/44 PEACEMAN-RACHFORD METHOD FOR NONLOCAL PROBLEMS 7 We denote Ix as an N ×N identity matrix and Iy as an (N − 1)× (N − 1) identity matrix. Now, we construct block matrices A1 =  Λx Λx . . . Λx Λx  , A2 =  2Ix −Ix . . . −Ix 2Ix −Ix . . . . . . . . . −Ix 2Ix −Ix . . . −Ix 2Ix  . (3.1) Every row and column of the matrices A1 and A2 has N−1 blocks of order N×N . We use Kronecker (tensor) product to investigate the eigenvalues and eigenvec- tors of the matrices A1 and A2. Definition 3.1 ([32]). Let A = {aij} and B = {bij} be two rectangular matrices of order m × n and p × q, accordingly. Matrix C of order mp × nq is called the Kronecker (tensor) product of matrices A and B C = A⊗B =  a11B a12B · · · a1nB a21B . . . . . . a2nB · · · . . . . . . · · · am1B am2B · · · anmB  . We rewrite the matrices A1 and A2 using Kronecker product A1 = Iy ⊗Λx, A2 = Λy ⊗ Ix. (3.2) Furthermore, we directly check that these matrices commute A1A2 = A2A1 = Λy ⊗Λx. (3.3) Now, we rewrite the systems of difference equations (2.13) and (2.14) with boundary conditions (2.11c) and (2.12c) in the matrix form( I + τ 2 A1 ) Un+1/2 = ( I− τ 2 A2 ) Un + τ 2 f̃n+1/2, (3.4)( I + τ 2 A2 ) Un+1 = ( I− τ 2 A1 ) Un+1/2 + τ 2 f̃n+1/2, (3.5) where I is (N − 1)N × (N − 1)N identity matrix and Un = { Unij } , f̃n+1/2 = {f̃n+1/2 ij }, i = 1, N, j = 1, N − 1. Expressing Un+1/2 from (3.4) and substituting into (3.5), we have Un+1 = SUn + τ 2 S1f̃ n+1/2, (3.6) where S = ( I + τ 2 A2 )−1 ( I− τ 2 A1 )( I + τ 2 A1 )−1 ( I− τ 2 A2 ) , (3.7) S1 = ( I + τ 2 A2 )−1 + ( I + τ 2 A2 )−1 ( I− τ 2 A1 )( I + τ 2 A2 )−1 . (3.8) 8 M. SAPAGOVAS, J. NOVICKIJ, R. ČIUPAILA EJDE-2022/44 We further define norms of any matrix and any vector to investigate stability of difference scheme (3.6). We also formulate theoretical proposition that we need to investigate stability of the difference method. Proposition 3.2 ([1, Th. 7.8]). Let ρ(A) be a spectral radius of an arbitrary square matrix A. If ε > 0 is given, then there exists a matrix norm ‖A‖∗ for which ‖A‖∗ ≤ ρ(A) + ε. Practically, we use a corollary of this proposition. Corollary 3.3. For any square matrix A, there exists a matrix norm ‖A‖∗ < 1 if and only if ρ(A) < 1. Now, we investigate when the condition ρ(S) < 1 (3.9) is valid for scheme (3.6). We find eigenvalues and eigenvectors of the matrices A1 and A2, as γ ∈ (0, 1). We formulate two eigenvalue problems (for the one- dimensional difference operator) ΛxV = µV, (3.10) where V = {Vi}, i = 1, N , and ΛyW = ηW, (3.11) where W = {Wj}, j = 1, N − 1. Suppose N is even. Then, the eigenvalues of the problem (3.10) are of the form (see [11]) µ1 = 4 h2 sin2 ψh 2 , µ2k = 4 h2 sin2 ( πk − ψh 2 ) , k = 1, N/2, µ2k+1 = 4 h2 sin2 ( πk + ψh 2 ) , k = 1, (N/2)− 1, (3.12) where ψ = arccos γ, 0 < ψ < 1. The corresponding eigenvectors V (k) := {V (k) i }, i = 1, N of the problem (3.10) are defined as V (1) i = sin(ψih), V (2k) i = sin ((2πk − ψ)ih) , k = 1, N/2, V (2k+1) i = sin ((2πk + ψ)ih) , k = 1, (N/2)− 1. (3.13) All the eigenvalues µ1, . . . , µN are positive and distinct and all the eigenvectors V (1), . . . , V (N) are linearly independent. If N is odd, expressions (3.12) and (3.13) are the same, except index k, which, in this case, is k = 1, (N − 1)/2. The eigenvalues and eigenvectors of difference problem (3.11) are of the form (see [24]) ηl = 4 h2 sin2 πlh 2 , l = 1, N − 1, (3.14) W (l) = { W (l) j } = {sin lπjh} , j, l = 1, N − 1. (3.15) All the eigenvalues (3.14) are positive and distinct and the eigenvectors (3.15) are linearly independent. EJDE-2022/44 PEACEMAN-RACHFORD METHOD FOR NONLOCAL PROBLEMS 9 Lemma 3.4. The matrices A1 and A2 have a common system of eigenvectors. Proof. Using eigenvectors V (k), k = 1, N of the matrix Λx and eigenvectors W (l), l = 1, N − 1 we construct new vector U (k,l) = W (l) ⊗ V (k) = ( W (l) 1 V (k),W (l) 2 V (k), . . . ,W (l) N−1V (k) )′ . (3.16) Using properties of tensor product (see [32]) we have (Iy ⊗Λx + Λy ⊗ Ix) (W ⊗V) = (µ+ η) (W ⊗V) . It means that the eigenvalues and the eigenvectors of a matrix A1 + A2 are µk + ηl and W (l) ⊗ V (k), respectively. Further, from the definitions of A1 and U (k,l) directly follows A1U (k,l) = µkU (k,l). (3.17) Using expressions (3.16) and (3.17) we have A2U (k,l) = ηlU (k,l). The eigenvectors W (l) ⊗ V (k) of the matrices A1 and A2 are the same. � Corollary 3.5 ([32]). Since the eigenvectors V (k), k = 1, N of the matrix Λx and the eigenvectors W (l), l = 1, N − 1 of the matrix Λy are linearly independent, then U (k,l) = W (l) ⊗ V (k) is the linearly independent system of vectors. Now, we define the norms of the matrices and vectors that we will use for inves- tigation of stability of scheme (3.6). We generate the matrix P, which columns are linearly independent eigenvectors U (k,l) of the matrix A1 (or A2). We define the norm of any m×m matrix A as ‖A‖∗ := ‖P−1AP‖2, (3.18) where ‖A‖2 = (max1≤i≤m λi(A ∗A))1/2 is the classical matrix norm and A∗ is the adjoint matrix. We define the compatible vector norm by the formula ‖U‖∗ = ‖P−1U‖2 = ( m∑ i=1 |Ũi|2 )1/2 , (3.19) where Ũi, i = 1,m are the coordinates of the vector P−1U. Indeed ‖AU‖∗ = ‖P−1AU‖2 = ‖P−1APP−1U‖2 ≤ ‖P−1AP‖2‖P−1U‖2 = ‖A‖∗‖U‖∗. For the special (symmetric or nonsymmetric) (m×m) matrix S we have relation ‖S‖∗ = ‖P−1SP‖2 = ‖J‖2 = max 1≤i≤m ‖µi(S)‖ = ρ(S), (3.20) where J = diag(µ1, . . . , µm), µi, i = 1,m are the eigenvalues and ρ(S) is the spectral radius of matrix S. One should not interpret (3.20) as the norm of any nonsymmetric matrix. This formula means that the norm of any matrix could be defined (see (3.18)) in a way that for particularly chosen matrix S equality ‖S‖∗ = ρ(S) is valid. We rewrite the vector norm (3.19) in an other form ‖U‖∗ = ( P−1U,P−1U )1/2 2 = ( (PP∗)−1U,U )1/2 2 = (DU,U) 1/2 2 , (3.21) where D is a positive definite matrix. The norm, defined as (3.21), is usually called energy norm, generated by the matrix D. 10 M. SAPAGOVAS, J. NOVICKIJ, R. ČIUPAILA EJDE-2022/44 Now, we return to the stability analysis of the scheme (3.6). We have ‖Un+1‖∗ = ‖S‖∗‖Un‖∗ + τ 2 ‖S1‖∗‖f̃n+1/2‖∗. Since the matrices A1 and A2 have the common system of eigenvectors, then, taking into account expressions (3.7) and (3.8), we have ‖S‖∗ = ρ(S) = max k,l ∣∣∣ (1− τ 2µk ) ( 1− τ 2ηl )( 1 + τ 2µk ) ( 1 + τ 2ηl ) ∣∣∣ := q < 1, where µk > 0 and ηl > 0. Similarly, ‖S1‖∗ = ρ(S1) = max k,l ∣∣∣ 1 1 + τ 2ηl + 1− τ 2µk( 1 + τ 2µk ) ( 1 + τ 2ηl ) ∣∣∣ < 1, where µk > 0 and ηl > 0. Therefore, we obtain the classical estimate for solution of the difference scheme ‖Un+1‖∗ ≤ q‖Un‖∗ + τ 2 ‖f̃n+1/2‖∗, 0 < q < 1, (3.22) which implies the stability of difference method in norm ‖U‖∗. So, the following theorem is valid. Theorem 3.6. If 0 < γ < 1, then the ADI method (2.13)–(2.14) for the differential problem (1.1)–(1.4) is stable in the norm ‖U‖∗. Proof. Proof of theorem follows from the above stated investigation. Really, if 0 < γ < 1, then all the eigenvalues µ1, . . . , µN are positive [11]. The eigenvalues η1, . . . , ηN are positive regardless of γ value. So, ρ(S) < 1 and ρ(S1) < 1. � Remark 3.7. In this article, the stability of ADI method, with nonsymmetric matrices A1 and A2, is based on the assumption that eigenvalues of these matrices are positive. To the authors’ knowledge, first time for the parabolic equations with nonlocal conditions this assumption to prove the stability of ADI method was used in [27]. However, several attempts to find suffiecient conditions for the stability or convergence of the ADI method with nonsymmetric matrices A1 and A2 were in the past (not for nonlocal boudary problems). One of these results is described in [24, Ch. 10, § 4.5]. It is proved that ‖S‖ < 1, if (Aαu, u) ≥ δα(u, u), δα > 0, α = 1, 2. It is easy to see that these assumptions are inappropriate for investigating the stability of ADI method presented in this article. All the eigenvalues of the matrix A1, defined as (3.1), are positive. Nevertheless, if u = {ui}, ui ≡ 1, then (A1u, u) = N∑ i,j=1 αij = 0. Therefore, assumption (Aαu, u) ≥ δα(u, u), δα > 0 is not fulfilled. 4. Numerical experiments Numerical experiments are performed to illustrate and confirm theoretical re- sults. We consider examples, where exact solution U of the problem (1.1)–(1.4) is not known. We compare two difference solutions U and U∗. U∗ is obtained twice EJDE-2022/44 PEACEMAN-RACHFORD METHOD FOR NONLOCAL PROBLEMS 11 reducing h and τ values. From theoretical investigation we know that approxima- tion error is O(h2 + τ2). So, it is presumed that if h and τ are reduced twice, maximum relative error ∆U := ∣∣∣max0≤i,j≤N U n ij −max0≤i,j≤N U ∗n ij max0≤i,j≤N Unij ∣∣∣, should be decreasing four times in each step. Case 1. γ = 0 (classical boundary condition), f(x, y, t) = 0, ϕ = y(1−y2)(π2 sin π 2x). h τ T = 1 T = 5 ∆U ∆U ratio ∆U ∆U ratio 0.1 0.1 1.07 · 10−6 1.25 · 10−6 0.05 0.05 2.81 · 10−7 3.81 3.29 · 10−7 3.8 0.025 0.025 7.26 · 10−8 3.87 8.55 · 10−8 3.85 0.0125 0.0125 1.85 · 10−8 3.92 2.20 · 10−8 3.89 0.00625 0.00625 4.66 · 10−9 3.97 5.58 · 10−9 3.94 Case 2. γ = 0.5, f(x, y, t) = 0, ϕ = y(1− y2)(π2 sin π 2x+ γ x 2 2 ). h τ T = 1 T = 5 ∆U ∆U ratio ∆U ∆U ratio 0.1 0.1 2.83 · 10−6 2.06 · 10−6 0.05 0.05 7.47 · 10−7 3.79 5.47 · 10−7 3.78 0.025 0.025 1.96 · 10−7 3.81 1.43 · 107 3.82 0.0125 0.0125 5.09 · 10−8 3.85 3.72 · 10−8 3.85 0.00625 0.00625 1.31 · 10−8 3.88 9.62 · 10−9 3.87 Case 3. γ = 1, f(x, y, t) = 0, ϕ = y(1− y2)(π2 sin π 2x+ γ x 2 2 ). h τ T = 1 T = 5 ∆U ∆U ratio ∆U ∆U ratio 0.1 0.1 3.54 · 10−6 2.88 · 10−6 0.05 0.05 1.01 · 10−6 3.51 8.14 · 10−7 3.54 0.025 0.025 2.80 · 10−7 3.60 2.27 · 10−7 3.59 0.0125 0.0125 7.51 · 10−8 3.73 6.19 · 10−8 3.66 0.00625 0.00625 1.93 · 10−8 3.89 1.64 · 10−8 3.78 Case 4. Unstable example γ = 2, f(x, y, t) = 0, ϕ = y(1− y2)(π2 sin π 2x+ γ x 2 2 ). h τ ∆U T = 1 T = 2 T = 5 0.025 0.025 1.80 · 10−4 1.80 · 10−3 3.01 · 10−2 0.025 0.0125 2.29 · 10−2 7.16 · 102 1.62 · 103 0.025 0.00625 3.17 · 102 4.76 · 104 6.12 · 104 0.025 0.003125 9.15 · 104 2.03 · 105 2.23 · 106 Remark 4.1. In this section we include examples with various γ values, starting from classical boundary conditions case (γ = 0) and some examples of nonlocal boundary conditions (γ = 0.5 and γ = 1). We also provide examples of experiments with different T values and show that stability does not depend on T choice. It follows from practical experiment that the constructed ADI method is stable and approximation error is of order O(h2 + τ2). In Case 4 the example of unstable scheme with γ = 2 is provided. It is known that if γ > 1, then the difference scheme for one-dimensional parabolic equation 12 M. SAPAGOVAS, J. NOVICKIJ, R. ČIUPAILA EJDE-2022/44 with nonlocal condition (1.2) is not stable with respect to the initial data [13]. It means that the second equation of the ADI method (2.13)–(2.14) is unstable. If the scheme is unstable, then ∆U grows while increasing the number of layers on the t axis. One can see from the experiment that either if τ decreases (h = const, T = const), or T increases (h = const, τ = const), then the relative error ∆U grows indefinitely. 5. Remarks and generalizations Stability of alternating direction method is proved in a special, quite complex vector norm ‖U‖∗ (see (3.19)). The equivalence of the same norm to the L2-norm for the one-dimensional parabolic equation with the nonlocal condition, similar to (1.2), is proved in [10]. To the authors’ knowledge, there was no investigation of stability, norms and spectrum for two-dimensional parabolic nonlocal problems of type (1.1)–(1.4). Nevertheless, introduced and investigated norm has important advantage. It follows from Proposition 3.2 and Corollary 3.3 that if difference method is not stable in norm ‖U‖∗, than ρ(S) > 1, and difference scheme is not stable in every other norm. This fact, for problems with nonlocal conditions, has been noticed and commented a lot previously (see, [13, 4, 29]). In this article, we theoretically investigated stability of alternating direction method in a case, when parameter γ in nonlocal condition is from interval (0, 1). With these γ values eigenvectors of both matrix S and matrix A1 +A2 are linearly independent, that is eigenvectors form a basis in vector space H(N−1)N . However, if γ = 1, then matrix Λx, the ADI method’s matrix S and Crank-Nicolson ma- trix A1 + A2 all have multiple eigenvalues corresponding to only one eigenvector. Therefore, these eigenvectors do not form basis in H(N−1)N (see [11]). This means that one cannot define norms ‖A‖∗ and ‖U‖∗ using formulas (3.18) and (3.19), because P−1 does not exist. Nevertheless, results of numerical experiments show that in sense of approximation accuracy variants γ ∈ (0, 1) and γ = 1 are both suitable. Not diving into the theoretical proofs, we explain this situation in two aspects. First, notice that Lemma 3.4 is valid for the case γ = 1 (see [32]). Therefore, according to Proposition 3.2 and Corollary 3.3, the same way as in case γ ∈ (0, 1), we have inequality ρ(S) < 1. From this inequality follows that norm ‖S‖∗ < 1 can be defined (without providing the method). In other words, in the case γ = 1 the ADI method is also stable. In this article we do not specify exact ‖S‖∗ and ‖U‖∗ expressions in the case γ = 1. This can be done at least in two ways. Just notice that the aim is not equality ‖S‖∗ = ρ(S), but ‖S‖∗ = ρ(S) + ε, ε < 0. One of the techniques is described in [25, Ch. 2, §3]. Another technique can be found in [15]. Second, in one-dimensional case problem (1.1)–(1.4) for both cases γ ∈ (0, 1) and γ = 1 has the same important property, related to concept of strong regularity of boundary conditions [20, 22]. Let γ = 1. Consider eigenvalue problem for nonlocal problem (1.1)–(1.4) in one-dimensional case. Define system of functions {ϕk(x)}, k = 1, 2, . . . for this eigenvalue problem, which consists of eigenvectors and associated vectors. Similar system {ψl(x)}, l = 1, 2, . . . consists of eigenvectors and associated vectors of adjoint eigenvalue problem. It is proved in [11] that both EJDE-2022/44 PEACEMAN-RACHFORD METHOD FOR NONLOCAL PROBLEMS 13 these systems are biorthonormal (ϕk, ψl) = ∫ 1 0 ϕk(x)ψl(x) dx = { 1, k = l, 0, k 6= l. In other words, the system of root functions of differential problem forms Risz basis in L2(0, 1), what is typical for problem with γ ∈ (0, 1) and, in general case, for problems with strongly regular boundary conditions. In more detail this property for two-dimensional hyperbolic equation with nonlocal condition is discussed in [22]. Furthermore, while investigating alternating direction method in this article, only stability is proved and nothing is said about convergence of difference method. We notice that similar situation occurs also in other papers, where difference meth- ods for differential equations are considered (see e.g. [11] and references therein). Often, as well as in this article, approximation of differential problem and stabil- ity of difference scheme with nonlocal conditions are considered in different vector norms. Additional investigation is required to prove equivalence of these norms. It is proved in one-dimensional case (see [10]) that the norm ‖U‖∗ is equivalent to vector L2-norm for all γ ∈ (0, 1) and γ = 1 values. This implies the convergence of the ADI method in L2-norm. To the authors’ knowledge, there are no investiga- tions for the two-dimensional case. It is obvious that matrices Λx, Λy, A1 and A2, defined in this article, are M -matrices [32]. Therefore, convergence of difference methods can be proved, by using properties of M -matrices [31, 6]. In the authors’ opinion, both of the ways to consider convergence of difference schemes with nonlocal conditions are worth separate investigation. References [1] K. E. Atkinson; An introduction to numerical analysis, John Wiley and Sons, 1989. [2] G. Berikelashvili, N. Khomeriki; On the convergence rate of a difference solution of the poisson equation with fully nonlocal constraints, Nonlinear Anal. Model. Control 19 (2014), no. 3, 367–381. [3] K. Bingelė, A. Bankauskienė, A. Štikonas; Spectrum curves for a discrete Sturm–Liouville problem with one integral boundary condition, Nonlinear Anal. Model. Control 24 (2019), no. 5, 755–774. [4] B. Cahlon, D. M. Kulkarni, P. Shi; Stepwise stability for the heat equation with a nonlocal constraint, SIAM J. Numer. Anal. 32 (1995), no. 2, 571–593. [5] J. R. Cannon The solution of the heat equation subject to specification of energy, Quart. Appl. Math. 21 (1963), no. 2, 155–160. [6] R. Čiupaila, M. Sapagovas, K. Pupalaigė; M-matrices and convergence of finite difference scheme for parabolic equation with an integral boundary condition, Math. Model. Anal. 25 (2020), no. 2, 167–183. [7] M. Dehghan; A new ADI technique for two-dimensional parabolic equation with an integral condition, Comput. Math. Applic. 43 (2002), 1477–1488. [8] M. Dehghan; Efficient techniques for the second-order parabolic equation subject to nonlocal specification, Applied Numer. Math. 52 (2005), no. 1, 39–62. [9] A. V. Gulin, N. I. Ionkin, V. A. Morozova; Difference schemes with nonlocal boundary con- ditions, Comput. Methods Appl. Math. 1 (2001), no. 1, 62–71. [10] A. V. Gulin, N. I. Ionkin, V. A. Morozova; Study of the norm in stability problems for nonlocal difference schemes, Differ. Equ. 42 (2006), no. 7, 974–984. [11] A. V. Gulin, N. I. Ionkin, V. A. Morozova; Ustojchivost nelokalnyh raznostnyh skhem (Sta- bility of nonlocal difference schemes), LKI, Moscow, 2008, (in Russian). [12] A. V. Gulin, V.A. Morozova; On the stability of a nonlocal finite-difference boundary value problem, Differ. Equ. 39 (2003), no. 7, 962–967. [13] A. V. Gulin, V. A. Morozova, N. S. Udovichenko; Stability criterion for a family of nonlocal difference schemes., Differ. Equ. 46 (2010), no. 7, 973–990. 14 M. SAPAGOVAS, J. NOVICKIJ, R. ČIUPAILA EJDE-2022/44 [14] N. I. Ionkin; The solution of a certain boundary value problem of the theory of heat conduction with a nonclassical boundary condition, Differ. Uravn. 13 (1977), no. 2, 294–304, (in Russian). [15] E. Isaacson, H. B. Keller; Analysis of numerical methods, Dover Publications, Inc., New York, 1994. [16] L. I. Kamynin; A boundary value problem in the theory of the heat conduction with nonclas- sical boundary condition, USSR Computational Mathematics and Mathematical Physics 4 (1964), no. 6, 33–59. [17] M. Kandemir, O. Sh. Mukhtarov; Nonlocal Sturm–Liouville problems with integral terms in the boundary conditions, Electr. J. Diff. Equat. 2017 (2017), no. 11, 1–12. [18] A. I. Kozhanov, L. S. Pul’kina; On the solvability of boundary value problems with a nonlocal boundary condition of integral form for multidimensional hyperbolic equations, Differ. Equ. 42 (2006), no. 9, 1233–1246. [19] J. Mart́ın-Vaquero, J. Vigo-Aguiar; On the numerical solution of the heat conduction equa- tions subject to nonlocal conditions, Applied Numer. Math. 59 (2009), no. 10, 2507–2514. [20] M. A. Naimark; Linejnye differencial’nye operatory (Linear differential operators), Moscow, Nauka, 1954, (in Russian). [21] A. M. Nakhushev; The equations of mathematical biology, Nauka, Moscow, 1995, (in Russian). [22] L. S. Pul’kina; Nonlocal problems for hyperbolic equations from the viewpoint of strongly regular boundary conditions, Electr. J. Diff. Equat. 2020 (2020), no. 28, 1–20. [23] B. P. Rynne; Spectral properties and nodal solutions for second-order, m-point, boundary value problems, Nonlinear Anal. Theory Methods Appl. 67 (2007), 3318–3327. [24] A. A. Samarskii; The theory of difference schemes, Marcel Dekker, Inc., New York, Basel, 2001. [25] A. A. Samarskii, A. V. Gulin; Chislennye metody (Numerical methods), Nauka, Moscow, 1989, (in Russian). [26] A. A. Samarskii, E. S. Nikolaev; Numerical methods for grid equations, Birkhüser, 1989. [27] M. Sapagovas, G. Kairytė, O. Štikoniené, A. Štikonas; Alternating direction method for a two-dimensional parabolic equation with a nonlocal boundary condition, Math. Model. Anal. 12 (2007), no. 1, 131–142. [28] M. Sapagovas, T. Meskauskas, and F. Ivanauskas; Numerical spectral analysis of a difference operator with non-local boundary conditions, Appl. Math. Comput. 218 (2012), no. 14, 7515– 7527. [29] M. Sapagovas, R. Čiupaila, Ž. Jokšienė, T. Meškauskas; Computational experiment for sta- bility analysis of difference schemes with nonlocal conditions, Informatica 24 (2013), no. 2, 275–290. [30] M. Sapagovas, A. Štikonas, O. Štikonienė; Alternating direction method for the Poisson equation with variable weight coefficients in an integral condition, Differ. Equ. 47 (2011), no. 8, 1176–1187. [31] M. Sapagovas, O. Štikonienė, K. Jakubelienė, R. Čiupaila; Finite difference method for bound- ary value problem for nonlinear elliptic equation with nonlocal conditions, Bound. Value Probl. 2019 (2019), no. 94, 1–16. [32] V. V. Voevodin, Y. A. Kuznecov; Matricy i vychisleniya (Matrices and computations), Nauka, Moscow, 1984, (in Russian). [33] A. Štikonas; A survey on stationary problems, Green’s functions and spectrum of Sturm– Liouville problem with nonlocal boundary conditions, Nonlinear Anal. Model. Control 19 (2014), no. 3, 301–334. [34] O. Štikonienė, M. Sapagovas, R. Čiupaila; On iterative methods for some elliptic equations with nonlocal conditions, Nonlinear Anal. Model. Control 19 (2014), no. 3, 517–535. Mifodijus Sapagovas Institute of Data Science and Digital Technologies, Vilnius University, Akademijos str. 4, LT-04812, Vilnius, Lithuania Email address: mifodijus.sapagovas@mif.vu.lt Jurij Novickij Institute of Data Science and Digital Technologies, Vilnius University, Akademijos str. 4, LT-04812 Vilnius, Lithuania Email address: jurij.novickij@mif.vu.lt EJDE-2022/44 PEACEMAN-RACHFORD METHOD FOR NONLOCAL PROBLEMS 15 Regimantas Čiupaila Vilnius Gediminas Technical University, Saulėtekio ave. 11, LT-10223, Vilnius, Lithua- nia Email address: regimantas.ciupaila@vilniustech.lt 1. Introduction 2. Numerical method and approximation error 3. Stability of difference scheme 4. Numerical experiments 5. Remarks and generalizations References