ADI Method for Pseudoparabolic Equation with Nonlocal Boundary Conditions

Abstract

This paper deals with the numerical solution of nonlocal boundary-value problem for two-dimensional pseudoparabolic equation which arise in many physical phenomena. A three-layer alternating direction implicit method is investigated for the solution of this problem. This method generalizes Peaceman–Rachford’s ADI method for the 2D parabolic equation. The stability of the proposed method is proved in the special norm. We investigate algebraic eigenvalue problem with nonsymmetric matrices to prove this stability. Numerical results are presented.

1. Introduction and Formulation of the Problem

In the past decades, the solution of the boundary-value problems with nonlocal boundary conditions has been an important and intensively investigated research area of numerical analysis and applied mathematics. We consider the third-order linear pseudoparabolic equation ($\eta >0$) with nonlocal integral conditions in the domain $D:=(0,L_x)\times (0,L_y)$

$$\frac{\partial u}{\partial t}=\frac{{\partial }^{2}u}{\partial x^2}+\frac{{\partial }^{2}u}{\partial y^2}+\eta \frac{\partial }{\partial t}\left(\frac{{\partial }^{2}u}{\partial x^2}+\frac{{\partial }^{2}u}{\partial y^2}\right)+f(x,y,t),(x,y,t)\in D\times (0,T],$$

(1)

subject to the initial and boundary conditions (BC)

$$\begin{array}{cc}\hfill u(0,y,t)& ={\gamma }_0{\int }_0^{L_x}u(x,y,t)\mathrm{d}x+v_l(y,t),(y,t)\in [0,L_y]\times [0,T],\hfill \end{array}$$

(2)

$$\begin{array}{cc}\hfill u(L_x,y,t)& ={\gamma }_1{\int }_0^{L_x}u(x,y,t)\mathrm{d}x+v_r(y,t),(y,t)\in [0,L_y]\times [0,T],\hfill \end{array}$$

(3)

$$\begin{array}{cc}\hfill u(x,0,t)& =v_b(x,t),u(x,L_y,t)=v_t(x,t),(x,t)\in [0,L_x]\times [0,T],\hfill \end{array}$$

(4)

$$\begin{array}{cc}\hfill u(x,y,0)& =u^0(x,y),(x,y)\in \overline{D}:=[0,L_x]\times [0,L_y].\hfill \end{array}$$

(5)

We have local Dirichlet BCs (4) in the y-direction and integral nonlocal boundary conditions (NBC) (2) and (3) in the x-direction.

Local and nonlocal boundary-value problems for the pseudoparabolic equations are the subject of intensive studies and a topic of great practical and theoretical interest, because many applied problems in physics, mechanics, and biology can be modelled using such equations.

For example, in the paper [1] the one-dimensional pseudoparabolic equation

$$\frac{\partial u}{\partial t}-\frac{{\partial }^{2}u}{\partial x^2}-\eta \frac{{\partial }^{3}u}{\partial x^2\partial t}=F(x,t,u),$$

(6)

$(x,t)\in (0,1)\times [0,T]$ subject to the initial conditions

$$u(x,0)=u_0\left(x\right),0\le x\le 1,$$

(7)

and to the integral conditions

$${\int }_0^{1}u(x,t)dx=E\left(t\right),{\int }_0^{1}xu(x,t)dx=G\left(t\right),0\le t\le T,$$

(8)

was considered. From a physical point of view, the problem (6)–(8) can be interpreted in the context of soil thermophysics. In this sense, (6) describes the dynamics of moisture in a subsoil layer $0<x<1$, while (8) represent the moisture moments [1,2].

Equations of type (6) with variable coefficients and additional terms also have many other applications in various physical situations, for example, in the theory of the two temperatures [3], in the study of the aggregation of population [4], or in the diffusion of imprisoned resonant radiation through a gas [5].

Theoretical research of pseudoparabolic equations with NBCs was started due to its applications for complex problems of science and technology. One of the first results for the third-order pseudoparabolic equations with NBCs for problems of soil dampness dynamics was considered in [6] (see also [1,7]). In the paper [8], the global existence of the weak solution for pseudoparabolic equation with the nonlocal source was considered.

A number of papers devoted to underground water flow dynamics modelled by pseudoparabolic equations with NBCs were published later, see [9,10]. Implicit finite difference schemes (FDS) for linear or nonlinear pseudoparabolic equations with Dirichlet BCs were introduced in the early 1970s [11,12]. Numerical methods for pseudoparabolic equations with NBCs have been of permanent interest for researchers during the last decades. Numerical methods for solving a one-dimensional nonlinear pseudoparabolic equation with integral conditions were introduced in [13,14]. FDS for linear 1D and 2D pseudoparabolic equations with various integral conditions were investigated in [15,16,17,18,19,20,21]. A separate class of pseudoparabolic equations consists of fractional pseudoparabolic equations, which were recently intensively studied [22,23,24].

Motivated by the works mentioned above, we study a class of two-dimensional pseudoparabolic equations with integral conditions (1)–(5). The goal of this paper is to generalize and investigate the Peaceman–Rachford alternating direction implicit (ADI) method [25] for pseudoparabolic Equation (1) with NBCs (2)–(5).

For the investigation of the stability of the ADI method, we study the structure of the spectrum of the corresponding difference operator. The structure of the spectrum of the differential and difference operators in the case of various nonlocal conditions were investigated in many papers (see, for example, [26,27,28,29]). Note that the presence of an integral term in the boundary condition can complicate the theoretical investigation of the numerical method. In this case, the matrix of the system of the difference equations is not symmetric or positive definite. The structure of the spectrum of such matrix is more complicated (and more substantive and informative) even in the case of Dirichlet or Neumann boundary value conditions. Furthermore, this structure strongly depends not only on the type of boundary conditions, but also on values of the parameters (or functions) included in the boundary conditions.

The first articles about the ADI method for the two-dimensional parabolic equations ($\eta =0$ in (1)) with NBCs showed that this method was also quite effective [30,31] in the case of NBC. In article [32], the proof of the stability of the ADI method for a parabolic equation was based on the analysis of the structure of the spectrum for a difference operator with nonlocal conditions. That approach was also relevant for the analysis of the convergence of the ADI method for elliptic equations with NBCs [33].

The presented ADI method is also suitable for solving an inverse problem for anomalous diffusion equation with a Riemann–Liouville derivative as well as parameter identification in fractional systems [34,35,36].

To the best of the authors’ knowledge, for the first time, a splitting method for pseudoparabolic Equation (1) with Dirichlet BCs was researched in articles [37,38].

The paper is organized as follows. In Section 2, we introduce the three-layer ADI method. Section 3 is devoted to the eigenvalue problem and in Section 4, we prove the main result of the present paper: the stability of the obtained three-layer difference scheme. In Section 5 the results of the numerical experiment are presented to demonstrate the accuracy and effectiveness of the finite difference scheme.

2. Alternating Direction Implicit (ADI) Method

Our goal is to construct a finite difference scheme for the pseudoparabolic problem and investigate its stability. The Peaceman–Rachford ADI method is used for developing a new ADI scheme for the third-order pseudoparabolic equation with NBC. The ADI method was introduced in 1955 by D.W. Peaceman and H.H. Rachford [25] and J. Douglas and H.H. Rachford [39] as a technique for the numerical solution of elliptic and parabolic differential equations. The theoretical and practical aspects of the ADI method led to extensions, generalizations, and ensuing applications far beyond the original application of a reservoir simulation. The advantage of the ADI method is that the equations that have to be solved at each step have a simpler structure and can be solved efficiently with the tridiagonal matrix algorithm. The ADI method is a predictor–corrector scheme where part of the difference operator is implicit in the initial prediction step and another part is implicit in the final correction step.

2.1. Notation

According to the standard technique of solving such problems by discretizations, we investigate finite difference schemes. We introduce grids with uniform steps ($1<N,M,L\in \mathbb{N}$):

$$\begin{array}{cc}\hfill {\overline{\omega }}_x^h& :=\{x_0=0,x_1,\cdots ,x_N=L_x\},h_x=x_i-x_{i-1}=L_x/N,i=\overline{1,N};\hfill \\ \hfill {\overline{\omega }}_y^h& :=\{y_0=0,y_1,\cdots ,y_M=L_y\},h_y=y_j-y_{j-1}=L_y/M,j=\overline{1,M};\hfill \\ \hfill {\overline{\omega }}^{\tau }& :=\{t_0=0,t_1,\cdots ,t_L=T\},\tau =t_k-t_{k-1}=T/L,k=\overline{1,L};\hfill \\ \hfill {\omega }_{1/2}^{\tau }& :=\{t_{1/2},\cdots ,t_{L-1/2}:t_{k-1/2}=(t_k+t_{k-1})/2,k=\overline{1,L}\};\hfill \\ \hfill {\omega }_x^h& :=\{x_1,\cdots ,x_{N-1}\};{\omega }_y^h:=\{y_1,\cdots ,y_{M-1}\};{\omega }^{\tau }:=\{t_1,\cdots ,t_L\}.\hfill \end{array}$$

We use the notation $U_{ij}^k:=U(x_i,y_j,t_k)$ for functions defined on the grid (or parts of this grid) $\overline{\omega }\times {\overline{\omega }}^{\tau }$ and $U_{ij}^{k-1/2}:=U(x_i,y_j,t_{k-1/2})$ on the grid $\overline{\omega }\times {\omega }_{1/2}^{\tau }$, where $\overline{\omega }:={\overline{\omega }}_x^h\times {\overline{\omega }}_y^h$. We omit indices if they are the same in the whole equation.

We use notation ${\mathbb{M}}_{k\times l}$ for the set of $k\times l$ matrices whose elements are real numbers, $k,l\in \mathbb{N}$. For any vector $\mathbf{u}\in {\mathbb{M}}_{m\times 1}$ and matrix $\mathbf{A}\in {\mathbb{M}}_{m\times m}$, discrete ${\ell }_2$-norms are defined as

$${\parallel \mathbf{u}\parallel }_2={\left(\sum _{i=1}^m{u}_i^2\right)}^{1/2},{\parallel \mathbf{A}\parallel }_2={\left(\underset{i}{max}\left|{\lambda }_i\left({\mathbf{A}}^{*}\mathbf{A}\right)\right|\right)}^{1/2},$$

(9)

where ${\mathbf{A}}^{*}$ is an adjoint matrix and ${\lambda }_i\left(\mathbf{A}\right)$ are eigenvalues of the matrix $\mathbf{A}$.

For functions $U_i=U\left(x_i\right)$ on the grid ${\overline{\omega }}_x^h$, we use the notation:

$$\begin{array}{cc}\hfill {\delta }_x^2{U}_i:=& \frac{U_{i-1}-2U_i+U_{i+1}}{h_x^2},(U,V):=\sum _{i=1}^{N-1}{U}_i{V}_i{h}_x\hfill \\ \hfill & [U,V]:=U_0{V}_0{h}_x/2+(U,V)+U_n{V}_n{h}_x/2.\hfill \end{array}$$

We introduce operators ${\delta }_x^2$, ${\delta }_y^2$ for a function on the grid $\overline{\omega }$:

$${\delta }_x^2{U}_{ij}:=\frac{U_{i-1,j}-2U_{ij}+U_{i+1,j}}{h_x^2},{\delta }_y^2{U}_{ij}:=\frac{U_{i,j-1}-2U_{ij}+U_{i,j+1}}{h_y^2},$$

and for the approximation integrals in the NBC, we define

$${[U,V]}_j:=U_{0j}{V}_{0j}{h}_x/2+{(U,V)}_j+U_{nj}{V}_{nj}{h}_x/2,{(U,V)}_j:=\sum _{i=1}^{N-1}{U}_{ij}{V}_{ij}{h}_x.$$

In a space H of grid functions $U_{ij}:=U(x_i,y_j)$ on the grid ${\omega }_x^h\times {\omega }_y^h$, we introduce an inner product

$${(U,V)}_H:=\sum _{j=1}^{M-1}{(U,V)}_j{h}_y.$$

Each such function is related to matrix $𝓤=\left(U_{ij}\right)\in {\mathbb{M}}_{(N-1)\times (M-1)}$. We choose one of the most obvious orderings and set a vector

$$\mathbf{U}:={(U_{11},\cdots ,U_{N-1,1},U_{12},\cdots ,U_{N-1,M-1})}^T=\mathrm{vec}\left(𝓤\right)\in {\mathbb{M}}_{(N-1)(M-1)\times 1}$$

for $U\in H$.

2.2. The ADI Method

Before writing down the solution method, we note that while using any splitting method for pseudoparabolic Equation (1), a one term with the third-order derivative $u_{txx}$ or $u_{tyy}$ has to be approximated in a lower layer, and for that we need at least two layers for each term. Therefore, we cannot write two-layer splitting FDS for two-dimensional pseudoparabolic equation, because the smallest number of layers is three.

Let us write the ADI method for pseudoparabolic Equation (1)

$$\begin{array}{ll}\hfill & \frac{U^{k+1/2}-U^k}{0.5\tau }={\delta }_x^2{U}^{k+1/2}+{\delta }_y^2{U}^k\hfill \\ \hfill & +\frac{\eta }{0.5\tau }\left({\delta }_x^2{U}^{k+1/2}-{\delta }_x^2{U}^k+{\delta }_y^2{U}^k-{\delta }_y^2{U}^{k-1/2}\right)+F^{k+1/2},\hfill \end{array}$$

(10)

$$\begin{array}{ll}\hfill & \frac{U^{k+1}-U^{k+1/2}}{0.5\tau }={\delta }_x^2{U}^{k+1/2}+{\delta }_y^2{U}^{k+1}\hfill \\ \hfill & +\frac{\eta }{0.5\tau }\left({\delta }_x^2{U}^{k+1/2}-{\delta }_x^2{U}^k+{\delta }_y^2{U}^{k+1}-{\delta }_y^2{U}^{k+1/2}\right)+F^{k+1/2}.\hfill \end{array}$$

(11)

In the case $\eta =0$, Equations (10) and (11) coincide with the classical Peaceman–Rachford method [25]. The truncation error of Equations (10) and (11) is $\mathcal{O}(\tau +h_x^2+h_y^2)$ if the solution of the differential problem (1)–(5) is smooth enough.

In order to start the calculation with this method, it is necessary to know the values of $U_{ij}^k$ in the two initial layers $k=0$ and $k=1/2$. The values of $U_{ij}^0$ are found from the initial condition (5). The values of $U_{ij}^{1/2}$ can be found by solving the system of difference equations written for the problem (1)–(5) using two layers. One can take ${\tilde{U}}_{ij}^{1/2}=U_{ij}^0$; in this case, we have error $|U_{ij}^{1/2}-{\tilde{U}}_{ij}^{1/2}|=\mathcal{O}\left(\tau \right)$.

Tridiagonal system (10) of order $N-1$ is solved for $j=\overline{1,M-1}$ with nonlocal boundary conditions

$$\begin{array}{cc}\hfill U^{k+1/2}{|}_{i=0}& ={\gamma }_0[1,U^{k+1/2}]+V_l^{k+1/2},\hfill \end{array}$$

(12)

$$\begin{array}{cc}\hfill U^{k+1/2}{|}_{i=N}& ={\gamma }_1[1,U^{k+1/2}]+V_r^{k+1/2}.\hfill \end{array}$$

(13)

The truncation error of the nonlocal boundary condition is $\mathcal{O}\left(h_x^2\right)$.

After the calculation of $U^{k+1/2}$, system (11) of order $M-1$ is solved for $i=\overline{1,N-1}$ with the Dirichlet BC

$$U^{k+1}{{|}_{j=0}=V_b^{k+1},U^{k+1}|}_{j=M}=V_t^{k+1}.$$

(14)

In order to investigate the stability of the ADI method (10)–(14), we rewrite this method in the most concise operator form. To this end, nonlocal conditions (12) and (13) are interpreted as a linear system with unknowns $U_{0j}^{k+1/2}$ and $U_{Nj}^{k+1/2}$ for every value of j. We express from conditions (12) and (13), the values $U_{0j}^{k+1/2}$ and $U_{Nj}^{k+1/2}$ via the remaining unknowns:

$$\begin{array}{cc}\hfill U^{k+1/2}{|}_{i=0}& ={\tilde{\gamma }}_0(1,U^{k+1/2})+{\tilde{V}}_l^{k+1/2},\hfill \end{array}$$

(15)

$$\begin{array}{cc}\hfill U^{k+1/2}{|}_{i=N}& ={\tilde{\gamma }}_1(1,U^{k+1/2})+{\tilde{V}}_r^{k+1/2},\hfill \end{array}$$

(16)

where ${\tilde{\gamma }}_0:={\gamma }_0{d}^{-1}$, ${\tilde{\gamma }}_1:={\gamma }_1{d}^{-1}$, ${\tilde{V}}_l:=(V_l+h_{x}c)d^{-1}$, ${\tilde{V}}_r:=(V_r-h_{x}c)d^{-1}$, $c:=({\gamma }_0{V}_r-{\gamma }_1{V}_l)/2$, $d:=1-\gamma h_x/2$, $\gamma :={\gamma }_0+{\gamma }_1$. We can write Formulas (15) and (16) in all cases when $d\ne 0$. Note, if $\gamma h_x<2$, then $d>0$.

Let us define square matrices ${\mathrm{\Lambda }}_x\in {\mathbb{M}}_{(N-1)\times (N-1)}$ and ${\mathrm{\Lambda }}_y\in {\mathbb{M}}_{(M-1)\times (M-1)}$ as

$$\begin{array}{cc}\hfill {\mathrm{\Lambda }}_x& :=\frac{1}{h_x^2}\left(\begin{array}{ccccc}2-{\tilde{\gamma }}_0{h}_x& -1-{\tilde{\gamma }}_0{h}_x& -{\tilde{\gamma }}_0{h}_x& \cdots & -{\tilde{\gamma }}_0{h}_x\\ -1& 2& -1& \cdots & 0\\ & \ddots & \ddots & \ddots & \\ & & -1& 2& -1\\ -{\tilde{\gamma }}_1{h}_x& \cdots & -{\tilde{\gamma }}_1{h}_x& -1-{\tilde{\gamma }}_1{h}_x& 2-{\tilde{\gamma }}_1{h}_x\end{array}\right),\hfill \\ \hfill {\mathrm{\Lambda }}_y& :=\frac{1}{h_y^2}\left(\begin{array}{ccccc}2& -1& & & \\ -1& 2& -1& & \\ & \ddots & \ddots & \ddots & \\ & & -1& 2& -1\\ & & & -1& 2\end{array}\right).\hfill \end{array}$$

Let ${\mathbf{I}}_x\in {\mathbb{M}}_{(N-1)\times (N-1)}$, ${\mathbf{I}}_y\in {\mathbb{M}}_{(M-1)\times (M-1)}$ be the identity matrices and $\mathbf{I}:={\mathbf{I}}_y\otimes {\mathbf{I}}_x$, where $\mathbf{A}\otimes \mathbf{B}$ denotes the Kronecker product of matrices $\mathbf{A}$ and $\mathbf{B}$. Then, we define matrices

$$\begin{array}{cc}\hfill {\mathbf{L}}_1& :={\mathbf{I}}_y\otimes {\mathrm{\Lambda }}_x=\mathrm{diag}({\mathrm{\Lambda }}_x,\cdots ,{\mathrm{\Lambda }}_x)\in {\mathbb{M}}_{(M-1)(N-1)\times (M-1)(N-1)},\hfill \\ \hfill {\mathbf{L}}_2& :={\mathrm{\Lambda }}_y\otimes {\mathbf{I}}_x=\left(\begin{array}{ccccc}2{\mathbf{I}}_y& -{\mathbf{I}}_y& & & \\ -{\mathbf{I}}_y& 2{\mathbf{I}}_y& -{\mathbf{I}}_y& & \\ & \ddots & \ddots & \ddots & \\ & & -{\mathbf{I}}_y& 2{\mathbf{I}}_y& -{\mathbf{I}}_y\\ & & & -{\mathbf{I}}_y& 2{\mathbf{I}}_y\end{array}\right)\in {\mathbb{M}}_{(M-1)(N-1)\times (M-1)(N-1)},\hfill \end{array}$$

which will be useful to rewrite the method (10) and (11), (14)–(16) in the operator form.

Lemma 1.

(see, [33]). The matrices ${\mathbf{L}}_1$ and ${\mathbf{L}}_2$ commute:

$${\mathbf{L}}_1{\mathbf{L}}_2={\mathbf{L}}_2{\mathbf{L}}_1={\mathrm{\Lambda }}_y\otimes {\mathrm{\Lambda }}_x.$$

(17)

Proof. 

From the properties of the Kronecker product, we have

$$\begin{array}{cc}\hfill {\mathbf{L}}_1{\mathbf{L}}_2& =({\mathbf{I}}_y\otimes {\mathrm{\Lambda }}_x)({\mathrm{\Lambda }}_y\otimes {\mathbf{I}}_x)=\left({\mathbf{I}}_y{\mathrm{\Lambda }}_y\right)\otimes \left({\mathrm{\Lambda }}_x{\mathbf{I}}_x\right)={\mathrm{\Lambda }}_y\otimes {\mathrm{\Lambda }}_x,\hfill \\ \hfill {\mathbf{L}}_1{\mathbf{L}}_2& =({\mathrm{\Lambda }}_y\otimes {\mathbf{I}}_x)({\mathbf{I}}_y\otimes {\mathrm{\Lambda }}_x)=\left({\mathrm{\Lambda }}_y{\mathbf{I}}_y\right)\otimes \left({\mathbf{I}}_x{\mathrm{\Lambda }}_x\right)={\mathrm{\Lambda }}_y\otimes {\mathrm{\Lambda }}_x.\hfill \end{array}$$

 □

Let us substitute expression (15) for $U_{0j}^{k+1/2}$ and expression (16) for $U_{nj}^{k+1/2}$ into Equation (10). Then, we can rewrite (10) in the matrix form

$${\mathbf{A}}_1{\mathbf{U}}^{k+1/2}+{\mathbf{B}}_1{\mathbf{U}}^k+{\mathbf{C}}_1{\mathbf{U}}^{k-1/2}=0.5\tau {\tilde{\mathbf{F}}}_1^{k+1/2},$$

(18)

where

$${\mathbf{A}}_1=\mathbf{I}+(\eta +0.5\tau ){\mathbf{L}}_1,{\mathbf{B}}_1=-\mathbf{I}+(\eta +0.5\tau ){\mathbf{L}}_2-\eta {\mathbf{L}}_1,{\mathbf{C}}_1=-\eta {\mathbf{L}}_2.$$

(19)

We rewrite Equation (11) using BC (14) in the form

$${\mathbf{A}}_2{\mathbf{U}}^{k+1}+{\mathbf{B}}_2{\mathbf{U}}^{k+1/2}+{\mathbf{C}}_2{\mathbf{U}}^k=0.5\tau {\tilde{\mathbf{F}}}_2^{k+1/2},$$

(20)

where

$${\mathbf{A}}_2=\mathbf{I}+(\eta +0.5\tau ){\mathbf{L}}_2,{\mathbf{B}}_2=-\mathbf{I}+(\eta +0.5\tau ){\mathbf{L}}_1-\eta {\mathbf{L}}_2,{\mathbf{C}}_2=-\eta {\mathbf{L}}_1.$$

(21)

We can write ${\tilde{\mathbf{F}}}_1^{k+1/2}$ and ${\tilde{\mathbf{F}}}_2^{k+1/2}$ in terms of known values ($F^{k+1/2}$, $V_b^{k+1}$, $V_t^{k+1}$, ${\tilde{V}}_l^{k+1/2}$, ${\tilde{V}}_r^{k+1/2}$), but the expressions of these functions are not important for the investigation of the stability of FDS. Each of Equations (18) and (20) separately corresponds to three-layer FDS (both together to four-layer FDS).

Lemma 2.

The matrices ${\mathbf{A}}_s$, ${\mathbf{B}}_s$, ${\mathbf{C}}_s$, and ${\mathbf{A}}_s^{-1}$, $s=1,2$, commute (if the inverse matrices exist).

Proof. 

The matrices ${\mathbf{L}}_1$ and ${\mathbf{L}}_2$ commute (see Lemma 1). Thus, the matrices ${\mathbf{A}}_s$, ${\mathbf{B}}_s$, and ${\mathbf{C}}_s$, defined by (19) and (21) also commute. If $\mathbf{A}\mathbf{B}=\mathbf{B}\mathbf{A}$ and ${\mathbf{A}}^{-\mathbf{1}}$ exists, then ${\mathbf{A}}^{-1}\mathbf{B}={\mathbf{A}}^{-1}\mathbf{B}\mathbf{A}{\mathbf{A}}^{-1}={\mathbf{A}}^{-1}\mathbf{A}\mathbf{B}{\mathbf{A}}^{-1}=\mathbf{B}{\mathbf{A}}^{-1}$. Therefore, it follows that matrices ${\mathbf{A}}_s$, ${\mathbf{B}}_s$, and ${\mathbf{C}}_s$ commute with ${\mathbf{A}}_s^{-1}$.  □

Remark 1.

We note that matrices ${\mathbf{A}}_1^{-1}$ and ${\mathbf{A}}_2^{-1}$ exist, as matrices ${\mathbf{A}}_1$ and ${\mathbf{A}}_2$ are strictly diagonally dominant.

Corollary 1.

The matrices ${\mathbf{A}}_1^{-1}{\mathbf{B}}_1$ and ${\mathbf{A}}_1^{-1}{\mathbf{C}}_1+\mu \mathbf{I}$ commute for all $\mu \in \mathbb{C}$.

Lemma 3.

The following equalities

$$\begin{array}{cc}\hfill & {\mathbf{A}}_s+{\mathbf{B}}_s+{\mathbf{C}}_s=0.5\tau ({\mathbf{L}}_1+{\mathbf{L}}_2)=2\mathbf{I}+{\mathbf{B}}_1+{\mathbf{B}}_2,s=1,2,\hfill \end{array}$$

(22)

$$\begin{array}{cc}\hfill & {\mathbf{A}}_s+{\mathbf{C}}_s=2\mathbf{I}+{\mathbf{B}}_{3-s},{\mathbf{A}}_s+{\mathbf{C}}_{3-s}=\mathbf{I}+0.5\tau {\mathbf{L}}_l,s=1,2,\hfill \end{array}$$

(23)

$$\begin{array}{cc}\hfill & ({\mathbf{A}}_1+{\mathbf{C}}_1)({\mathbf{A}}_2+{\mathbf{C}}_2)={\mathbf{B}}_1{\mathbf{B}}_2+\tau ({\mathbf{L}}_1+{\mathbf{L}}_2)\hfill \end{array}$$

(24)

are valid.

Proof. 

The equalities (22) and (23) are obvious. Finally,

$$\begin{array}{cc}\hfill ({\mathbf{A}}_1+{\mathbf{C}}_1)({\mathbf{A}}_2+{\mathbf{C}}_2)& =(2\mathbf{I}+{\mathbf{B}}_2)(2\mathbf{I}+{\mathbf{B}}_1)=4\mathbf{I}+2{\mathbf{B}}_1+2{\mathbf{B}}_2+{\mathbf{B}}_2{\mathbf{B}}_1\hfill \\ \hfill & ={\mathbf{B}}_2{\mathbf{B}}_1+\tau ({\mathbf{L}}_1+{\mathbf{L}}_2)={\mathbf{B}}_1{\mathbf{B}}_2+\tau ({\mathbf{L}}_1+{\mathbf{L}}_2).\hfill \end{array}$$

 □

Corollary 2.

The following equality

$$-{\mathbf{B}}_1{\mathbf{B}}_2+{\mathbf{A}}_1{\mathbf{C}}_2+{\mathbf{A}}_2{\mathbf{C}}_1=-{\mathbf{A}}_1{\mathbf{A}}_2-{\mathbf{C}}_1{\mathbf{C}}_2+\tau ({\mathbf{L}}_1+{\mathbf{L}}_2)$$

is valid.

The methodology of the investigation of the stability of three-layer schemes for the second-order parabolic equation with Dirichlet boundary conditions is created in monograph [40]. According to that methodology, the three-layer difference scheme (18) or (20) must be rewritten in the canonical form. For example, the canonical form for (18) is

$$\tilde{\mathbf{B}}\frac{{\mathbf{U}}^{k+1/2}-{\mathbf{U}}^{k-1/2}}{\tau }+\tilde{\mathbf{R}}({\mathbf{U}}^{k+1/2}-2{\mathbf{U}}^k+{\mathbf{U}}^{k-1/2})+\tilde{\mathbf{A}}{\mathbf{U}}^k={\tilde{\mathbf{F}}}_1^{k+1/2},$$

where

$$\tilde{\mathbf{B}}=\mathbf{I}+(\eta +0.5\tau ){\mathbf{L}}_1+\eta {\mathbf{L}}_2,\tau \tilde{\mathbf{R}}=\mathbf{I}+(\eta +0.5\tau ){\mathbf{L}}_1-\eta {\mathbf{L}}_2,\tilde{\mathbf{A}}={\mathbf{L}}_1+{\mathbf{L}}_2.$$

For the stability of difference scheme (18), the matrices $\tilde{\mathbf{R}}$ and $\tilde{\mathbf{A}}$ must be positive definite [40]. However, in the case of nonlocal conditions (3) and (4), these matrices are nonsymmetric because of nonsymmetrical matrix ${\mathrm{\Lambda }}_x$. Furthermore, the symmetric matrix $\tilde{\mathbf{R}}$ is not positive definite in the case of pseudoparabolic equation ($\eta \ne 0$). Thus, we investigate the stability of difference scheme (18) by using the other approach.

2.3. Reduction of the Three-Layer Scheme to a Two-Layer Scheme

In order to reduce the three-layer difference scheme to a two-layer system, we define vectors in ${\mathbb{M}}_{2(N-1)(M-1)\times 1}$:

$${\mathbf{Y}}^k=\left(\begin{array}{c}{\mathbf{U}}^k\\ {\mathbf{U}}^{k-1/2}\end{array}\right),{\mathbf{F}}_s^{k+1/2}=\left(\begin{array}{c}\frac{\tau }{2}{\mathbf{A}}_s^{-1}{\tilde{\mathbf{F}}}_s^{k+1/2}\\ \mathbf{0}\end{array}\right),s=1,2.$$

Then, we rewrite Equations (18) and (20) as

$$\begin{array}{c}\hfill {\mathbf{Y}}^{k+1/2}={\mathbf{S}}_1{\mathbf{Y}}^k+{\mathbf{F}}_1^{k+1/2},{\mathbf{Y}}^{k+1}={\mathbf{S}}_2{\mathbf{Y}}^{k+1/2}+{\mathbf{F}}_2^{k+1/2},\end{array}$$

where

$${\mathbf{S}}_l=\left(\begin{array}{cc}-{\mathbf{A}}_l^{-1}{\mathbf{B}}_l& -{\mathbf{A}}_l^{-1}{\mathbf{C}}_l\\ \mathbf{I}& \mathbf{0}\end{array}\right)\in {\mathbb{M}}_{2(N-1)(M-1)\times 2(N-1)(M-1)},l=1,2.$$

From this, we obtain

$${\mathbf{Y}}^{k+1}=\mathbf{S}{\mathbf{Y}}^k+{\mathbf{F}}^{k+1/2},$$

(25)

where

$$\mathbf{S}={\mathbf{S}}_2{\mathbf{S}}_1=\left(\begin{array}{cc}{\mathbf{A}}_2^{-1}{\mathbf{B}}_2{\mathbf{A}}_1^{-1}{\mathbf{B}}_1-{\mathbf{A}}_2^{-1}{\mathbf{C}}_2& {\mathbf{A}}_2^{-1}{\mathbf{B}}_2{\mathbf{A}}_1^{-1}{\mathbf{C}}_1\\ -{\mathbf{A}}_1^{-1}{\mathbf{B}}_1& -{\mathbf{A}}_1^{-1}{\mathbf{C}}_1\end{array}\right),$$

(26)

and ${\mathbf{F}}^{k+1/2}={\mathbf{S}}_2{\mathbf{F}}_1^{k+1/2}+{\mathbf{F}}_2^{k+1/2}$.

3. Eigenvalues of Matrix $\mathbf{S}$

3.1. Eigenvalues of the Matrices ${\mathrm{\Lambda }}_x$ and ${\mathrm{\Lambda }}_y$

It is known that the eigenvalues of ${\mathrm{\Lambda }}_y$ are positive and the corresponding eigenvectors are orthogonal and linearly independent:

$${\lambda }_m\left({\mathrm{\Lambda }}_y\right)=\frac{4}{h_y^2}{sin}^2\frac{\pi m{h}_y}{2L_y},{\mathbf{W}}^m={(W_1^m,\cdots ,W_{M-1}^m)}^T,m=\overline{1,M-1},$$

where

$$W_j^m=sin\frac{\pi mj{h}_y}{L_y},j=\overline{1,M-1}.$$

The eigenvalue problem for matrix ${\mathrm{\Lambda }}_x$ is equivalent to problem

$$\begin{array}{cc}\hfill -{\delta }_x^2{V}_i=\lambda V_i,& i=\overline{1,N-1},\hfill \\ \hfill V_0={\gamma }_0[1,V],& V_n={\gamma }_1[1,V].\hfill \end{array}$$

All eigenvalues ${\lambda }_n\left({\mathrm{\Lambda }}_y\right)$, $n=\overline{1,N-1}$, of this problem are real numbers and eigenvectors

$${\mathbf{V}}^n={(V_1^n,\cdots ,V_{N-1}^n)}^T,n=\overline{1,N-1},$$

are linearly independent. Furthermore, if ${\gamma }_0+{\gamma }_1<2$, then all eigenvalues ${\lambda }_n\left({\mathrm{\Lambda }}_x\right)$, $n=\overline{1,N-1}$, are positive, more precisely, ${\lambda }_n\left({\mathrm{\Lambda }}_x\right)\in (0,4/h_x^2)$ [28].

3.2. Eigenvalues and Eigenvectors of the Matrices ${\mathbf{L}}_1$ and ${\mathbf{L}}_2$

Vectors

$${\mathbf{U}}^{m,n}={\mathbf{W}}^m\otimes {\mathbf{V}}^n,n=\overline{1,N-1},m=\overline{1,M-1},$$

(27)

are linearly independent in the vector space ${\mathbb{R}}^{2(N-1)(M-1)}$ according to the Kronecker product properties [41].

Since ${\mathbf{L}}_1={\mathbf{I}}_y\otimes {\mathrm{\Lambda }}_x$ and ${\mathbf{L}}_2={\mathrm{\Lambda }}_y\otimes {\mathbf{I}}_x$, we have

$$\begin{array}{cc}\hfill ({\mathbf{I}}_y\otimes {\mathrm{\Lambda }}_x)({\mathbf{W}}^m\otimes {\mathbf{V}}^n)& =\left({\mathbf{I}}_y{\mathbf{W}}^m\right)\otimes \left({\mathrm{\Lambda }}_x{\mathbf{V}}^n\right)={\lambda }_n\left({\mathrm{\Lambda }}_x\right){\mathbf{W}}^m\otimes {\mathbf{V}}^n,\hfill \\ \hfill ({\mathrm{\Lambda }}_y\otimes {\mathbf{I}}_x)({\mathbf{W}}^m\otimes {\mathbf{V}}^n)& =\left({\mathrm{\Lambda }}_y{\mathbf{W}}^m\right)\otimes \left({\mathbf{I}}_x{\mathbf{V}}^n\right)={\lambda }_m\left({\mathrm{\Lambda }}_y\right){\mathbf{W}}^m\otimes {\mathbf{V}}^n,\hfill \end{array}$$

and we get

$${\mathbf{L}}_1{\mathbf{U}}^{m,n}={\lambda }_n\left({\mathrm{\Lambda }}_x\right){\mathbf{U}}^{m,n},{\mathbf{L}}_2{\mathbf{U}}^{m,n}={\lambda }_m\left({\mathrm{\Lambda }}_y\right){\mathbf{U}}^{m,n}.$$

(28)

Therefore, matrices ${\mathbf{L}}_1$ and ${\mathbf{L}}_2$ have the same system of eigenvectors. The eigenvalues of matrix ${\mathbf{L}}_1$ are ${\lambda }_n^{\left(1\right)}={\lambda }_n\left({\mathrm{\Lambda }}_x\right)$, $n=\overline{1,N-1}$, and the geometric multiplicity of each eigenvalue is $M-1$; the eigenvalues of matrix ${\mathbf{L}}_2$ are ${\lambda }_m^{\left(2\right)}={\lambda }_m\left({\mathrm{\Lambda }}_y\right)$, $m=\overline{1,M-1}$, and the geometric multiplicity of each eigenvalue is $N-1$.

Lemma 4.

The matrices ${\mathbf{A}}_s$, ${\mathbf{B}}_s$, ${\mathbf{C}}_s$, and ${\mathbf{L}}_s$, $s=1,2,$ commute and have a common system of linearly independent eigenvectors.

Proof. 

The matrices ${\mathrm{\Lambda }}_1$ and ${\mathrm{\Lambda }}_2$ commute and have the same system of eigenvectors. Thus, according to (19) and (21), vectors ${\mathbf{U}}^{m,n}$ are eigenvectors of matrices ${\mathbf{A}}_s$, ${\mathbf{B}}_s$, ${\mathbf{C}}_s$, too.  □

3.3. Eigenvalues of the Matrix $\mathbf{S}$

The equation for the eigenvalues of matrix $\mathbf{S}$ is

$$\left|\begin{array}{cc}{\mathbf{A}}_2^{-1}{\mathbf{B}}_2{\mathbf{A}}_1^{-1}{\mathbf{B}}_1-{\mathbf{A}}_2^{-1}{\mathbf{C}}_2-\mu \mathbf{I}& {\mathbf{A}}_2^{-1}{\mathbf{B}}_2{\mathbf{A}}_1^{-1}{\mathbf{C}}_1\\ -{\mathbf{A}}_1^{-1}{\mathbf{B}}_1& -{\mathbf{A}}_1^{-1}{\mathbf{C}}_1-\mu \mathbf{I}\end{array}\right|=0.$$

The determinant on the left side is equal to

$$\left|\begin{array}{cc}-{\mathbf{A}}_2^{-1}{\mathbf{C}}_2-\mu \mathbf{I}& -\mu {\mathbf{A}}_2^{-1}{\mathbf{B}}_2\\ -{\mathbf{A}}_1^{-1}{\mathbf{B}}_1& -{\mathbf{A}}_1^{-1}{\mathbf{C}}_1-\mu \mathbf{I}\end{array}\right|=\frac{1}{det\left({\mathbf{A}}_1{\mathbf{A}}_2\right)}\left|\begin{array}{cc}{\mathbf{C}}_2+\mu {\mathbf{A}}_2& \mu {\mathbf{B}}_2\\ {\mathbf{B}}_1& {\mathbf{C}}_1+\mu {\mathbf{A}}_1\end{array}\right|.$$

Using Lemma 2 and Corollary 1, we get equation

$$det\left({\mu }^2\mathbf{A}+\mu \mathbf{B}+\mathbf{C}\right)=0,$$

(29)

where

$$\begin{array}{cc}\hfill \mathbf{A}& ={\mathbf{A}}_1{\mathbf{A}}_2,\mathbf{C}={\mathbf{C}}_1{\mathbf{C}}_2,\hfill \end{array}$$

(30)

$$\begin{array}{cc}\hfill \mathbf{B}& =-{\mathbf{B}}_1{\mathbf{B}}_2+{\mathbf{A}}_1{\mathbf{C}}_2+{\mathbf{A}}_2{\mathbf{C}}_1=-{\mathbf{A}}_1{\mathbf{A}}_2-{\mathbf{C}}_1{\mathbf{C}}_2+\tau ({\mathbf{L}}_1+{\mathbf{L}}_2).\hfill \end{array}$$

(31)

The last equality follows from Corollary 2. Equation (29) is the characteristic equation of the general nonlinear eigenvalue problem

$${\mu }^2\mathbf{A}\mathbf{U}+\mu \mathbf{B}\mathbf{U}+\mathbf{C}\mathbf{U}=0.$$

(32)

Thus, the second order eigenvalue problem has $2(N-1)(M-1)$ eigenvalues. As a result, the next lemma is valid.

Lemma 5.

The eigenvalues of matrix $\mathbf{S}$ coincide with the eigenvalues of the nonlinear eigenvalue problem (32).

The nonlinear eigenvalue problem of such type when matrices $\mathbf{A}$, $\mathbf{B}$, and $\mathbf{C}$ are symmetric is well known and has been considered in many works (see, [42]). In our case, these matrices are nonsymmetric, but it is possible to investigate nonlinear eigenvalue problem (32) using another useful property of these matrices. Namely, matrices $\mathbf{A}$, $\mathbf{B}$, and $\mathbf{C}$ have the same system of eigenvectors.

We can find the eigenvalues of nonlinear problem (32). Substituting $\mathbf{U}={\mathbf{U}}^{m,n}$ into Equation (32) and taking into account that ${\mathbf{U}}^{m,n}\not\equiv \mathbf{0}$, we obtain

$${\mu }^2{\lambda }_{mn}\left(\mathbf{A}\right)+\mu {\lambda }_{mn}\left(\mathbf{B}\right)+{\lambda }_{mn}\left(\mathbf{C}\right)=0,m=\overline{1,M-1},n=\overline{1,N-1}.$$

(33)

The eigenvalues of matrix $\mathbf{S}$ are roots of Equation (33). We rewrite (33) as

$$a_{mn}{\mu }^2+b_{mn}\mu +c_{mn}=0,m=\overline{1,M-1},n=\overline{1,N-1},$$

(34)

where

$$\begin{array}{cc}\hfill a_{mn}& =\left(1+(\eta +0.5\tau ){\lambda }_n^{\left(1\right)}\right)\left(1+(\eta +0.5\tau ){\lambda }_m^{\left(2\right)}\right),\hfill \end{array}$$

(35)

$$\begin{array}{cc}\hfill b_{mn}& =-a_{mn}-c_{mn}+\tau \left({\lambda }_n^{\left(1\right)}+{\lambda }_m^{\left(2\right)}\right),\hfill \end{array}$$

(36)

$$\begin{array}{cc}\hfill c_{mn}& ={\eta }^2{\lambda }_n^{\left(1\right)}{\lambda }_m^{\left(2\right)}.\hfill \end{array}$$

(37)

If ${\lambda }_n^{\left(1\right)}>0$, ${\lambda }_m^{\left(2\right)}>0$, then taking into account $\eta >0$, we get that $c_{mn}>0$. From (35), we get

$$a_{mn}>c_{mn},a_{mn}>0.5\tau \left({\lambda }_n^{\left(1\right)}+{\lambda }_m^{\left(2\right)}\right).$$

Thus, we estimate

$$b_{mn}<-a_{mn}-c_{mn}+2a_{mn}=a_{nm}-c_{mn}<a_{nm}+c_{mn}.$$

From (36), we have $-b_{mn}<a_{mn}+c_{mn}$. Finally, we prove

$$|c_{mn}|<a_{mn},\left|b_{mn}\right|<c_{mn}+a_{mn}.$$

Lemma 6.

If ${\gamma }_0+{\gamma }_1<2$, then for the roots of Equation (34), we have $\left|\mu \right|<1$.

Proof. 

If ${\gamma }_0+{\gamma }_1<2$, then ${\lambda }_n\left({\mathrm{\Lambda }}_x\right)>0$, $n=\overline{1,N-1}$. Eigenvalues ${\lambda }_m\left({\mathrm{\Lambda }}_y\right)$, $m=\overline{1,M-1}$, are positive. Thus, we have ${\lambda }_n^{\left(1\right)}>0$, ${\lambda }_m^{\left(2\right)}>0$ for $m=\overline{1,M-1}$, $n=\overline{1,N-1}$. According to Hurwitz’s criterion [43], the roots of polynomial ${\mu }^2+b\mu +c$ satisfy condition $\left|\mu \right|<1$ if and only if $\left|c\right|<1$, $\left|b\right|<c+1$. Therefore, the conditions of Hurwitz’s criterion are fulfilled for $b=a_{mn}^{-1}{b}_{mn}$, $c=a_{mn}^{-1}{c}_{mn}$, and we get that $\left|\mu \right|<1$.  □

Remark 2.

In the case $\eta =0$, we have $c_{mn}=0$, and the roots are ${\mu }_1=0$, ${\mu }_2=-b_{mn}{a}_{mn}^{-1}=\left(1-0.5\tau {\lambda }_n^{\left(1\right)}\right)\left(1-0.5\tau {\lambda }_m^{\left(2\right)}\right){\left(1+0.5\tau {\lambda }_n^{\left(1\right)}\right)}^{-1}{\left(1+0.5\tau {\lambda }_m^{\left(2\right)}\right)}^{-1}$, and $|{\mu }_2|<1$.

In the limit case $\tau =0$, we have equation $a_{mn}{\mu }^2-(a_{mn}+c_{mn})\mu +c_{mn}=0$. The roots of this equation are ${\mu }_1=c_{mn}{a}_{mn}^{-1}={\eta }^2{\lambda }_n^{\left(1\right)}{\lambda }_m^{\left(2\right)}{\left(1+\eta {\lambda }_n^{\left(1\right)}\right)}^{-1}{\left(1+\eta {\lambda }_m^{\left(2\right)}\right)}^{-1}$, ${\mu }_2=1$, and $|{\mu }_1|<1$.

4. Investigation of the Stability of Finite Difference Scheme

We examine, based on Lemma 6, the stability of the ADI method (10)–(11) or otherwise the difference scheme (25).

From Lemma 6, we can formulate the following proposition.

Proposition 1.

If $\left|\mu \right|<1$, where μ is the eigenvalue of matrix $\mathbf{S}$, then the difference scheme ${\mathbf{Y}}^{n+1}=\mathbf{S}{\mathbf{Y}}^n+{\mathbf{F}}^{n+1/2}$ is stable regardless of whether the matrix $\mathbf{S}$ is symmetric or not.

This proposition on the stability of the two-layer difference scheme for a parabolic equation has been known for a long time and was formulated before solving problems with nonlocal conditions (see, [44]). In order to understand the influence of nonlocal conditions on the stability of the difference scheme, we study in detail the case when $\mathbf{S}$ is a nonsymmetric matrix.

If ${\gamma }_1={\gamma }_2=0$ in problem (10)–(13), then matrix $\mathbf{S}$ is symmetric. In this case, all eigenvalues of $\mathbf{S}$ are real and it follows from $\left|\mu \right|<1$ that

$$\varrho \left(\mathbf{S}\right):=max\left|\mu \right(\mathbf{S}\left)\right|<1.$$

(38)

If $\mathbf{S}$ is symmetric, we have

$${\parallel \mathbf{S}\parallel }_2={\left(\underset{i}{max}\left|{\lambda }_i\left({\mathbf{S}}^2\right)\right|\right)}^{1/2}=\varrho \left(\mathbf{S}\right).$$

(39)

Thus, the condition of the stability of the finite difference scheme (25) can be formulated as the following theorem.

Theorem 1.

If ${\gamma }_1={\gamma }_2=0$, then

$${\parallel \mathbf{S}\parallel }_2<1$$

(40)

and ADI method (9)–(12) is stable in the discrete vector ${\ell }_2$-norm.

Let us assume on the contrary that either ${\gamma }_1$ or ${\gamma }_2$ is nonzero in problem (1)–(5). Thus, we investigate the stability of difference scheme (10)–(13) or the system (25) with nonsymmetric matrix $\mathbf{S}$. In this case, the equality (39) does not follow from the definition of the norm (9). That is, for nonsymmetric matrix $\mathbf{S}$, the spectral radius $\varrho \left(\mathbf{S}\right)$ is not the matrix norm. For this, we use the statement from linear algebra about the norm of the matrix.

Proposition 2.

(see, [45], Th. 7.8). Let $\varrho \left(\mathbf{A}\right)$ be a spectral radius of an arbitrary square matrix $\mathbf{A}$. If $\epsilon >0$ is given, then there exists a matrix norm ${\parallel \mathbf{A}\parallel }_{*}$ for which

$${\parallel \mathbf{A}\parallel }_{*}\le \varrho \left(\mathbf{A}\right)+\epsilon .$$

(41)

We can formulate some corollary of this proposition.

Corollary 3.

For any square matrix $\mathbf{A}$, there exists a matrix norm ${\parallel \mathbf{A}\parallel }_{*}<1$ if and only if $\varrho \left(\mathbf{A}\right)<1$.

However, the proof of Proposition 2, as well as a construction of norm ${\parallel \mathbf{A}\parallel }_{*}$, is nontrivial (see, [46] p. 12, or [47] Chapter 11.2, Section 3.4).

Now, we assume that matrix $\mathbf{S}$ has a system of linearly independent eigenvectors. In this case, we use the following proposition from linear algebra [48].

Proposition 3.

Let $\parallel \mathbf{A}\parallel$ and $\parallel \mathbf{u}\parallel$ be compatible matrix and vector norms and $\mathbf{P}$ be a nonsingular matrix ($det\mathbf{P}\ne 0$). Then,

$${\parallel \mathbf{A}\parallel }_{*}=\parallel {\mathbf{P}}^{-1}\mathbf{A}\mathbf{P}\parallel$$

(42)

and

$${\parallel \mathbf{u}\parallel }_{*}=\parallel {\mathbf{P}}^{-1}\mathbf{u}\parallel$$

(43)

are also compatible matrix and vector norms.

Using the assumption that the eigenvectors of matrix $\mathbf{S}$ are linearly independent, we define matrix $\mathbf{P}$, whose columns are eigenvectors of $\mathbf{S}$. We use the ${\ell }_2$-norm ${\parallel ·\parallel }_2$ in (42) and (43). Therefore, we get

$${\parallel \mathbf{u}\parallel }_{*}={\parallel {\mathbf{P}}^{-1}\mathbf{u}\parallel }_2={({\mathbf{P}}^{-1}\mathbf{u},{\mathbf{P}}^{-1}\mathbf{u})}^{1/2}={(\mathbf{D}\mathbf{u},\mathbf{u})}^{1/2},\mathbf{D}={\left(\mathbf{P}{\mathbf{P}}^{*}\right)}^{-1},$$

(44)

$${\parallel \mathbf{S}\parallel }_{*}=\parallel {\mathbf{P}}^{-1}{\mathbf{S}\mathbf{P}\parallel }_2={\parallel \mathbf{J}\parallel }_2=\underset{i}{max}\left|{\mu }_i\left(\mathbf{S}\right)\right|=\varrho \left(\mathbf{S}\right),$$

(45)

where $\mathbf{D}$ is a positive defined matrix and $\mathbf{J}$ is the Jordan form of $\mathbf{S}$.

Now, we can generalize Theorem 1.

Theorem 2.

If ${\gamma }_1+{\gamma }_2<2$, then the ADI method (10) and (11) with nonlocal conditions (12) and (13) is stable in some vector norm. If additionally, the eigenvectors of matrix $\mathbf{S}$ are linearly independent, then the ADI method is stable in norm ${\parallel \mathbf{u}\parallel }_{*}={(\mathbf{D}\mathbf{u},\mathbf{u})}^{1/2}$, generated by the self-adjoint positive defined operator (matrix) $\mathbf{D}={\left(\mathbf{P}{\mathbf{P}}^{*}\right)}^{-1}$, where $\mathbf{P}$ is a matrix whose columns are linearly independent eigenvectors of $\mathbf{S}$.

Remark 3.

As far as the authors know, the eigenvectors of a matrix $\mathbf{S}$ are linearly independent if there are some additional conditions. For example, in the case when the eigenvalues of $\mathbf{S}$ are real. In this article, we leave this question open.

5. Numerical Examples

One of the aims of our numerical simulations was to demonstrate the theoretical results obtained in the previous sections. Two numerical examples illustrate the effectiveness of the present ADI scheme. Another aim was to investigate numerically the influence of parameter $\tau$, T, and $\eta$ on the accuracy of the solution.

Problem 1.

We consider the model problem in the domain $D:=(0,1)\times (0,1)$

$$\frac{\partial u}{\partial t}=\frac{{\partial }^{2}u}{\partial x^2}+\frac{{\partial }^{2}u}{\partial y^2}+\eta \frac{\partial }{\partial t}\left(\frac{{\partial }^{2}u}{\partial x^2}+\frac{{\partial }^{2}u}{\partial y^2}\right)+f(x,y,t),(x,y,t)\in D\times (0,T],$$

where

$$f(x,y,t)={\mathrm{e}}^{-t}sin\left(\pi x\right)sin\left(\pi y\right)\left(-1+2{\pi }^2(1-\eta )\right),$$

subject to the initial and boundary conditions

$$\begin{array}{cc}\hfill u(0,y,t)& ={\gamma }_0{\int }_0^{1}u(x,y,t)\mathrm{d}x-2\frac{{\gamma }_0}{\pi }{\mathrm{e}}^{-t}sin\left(\pi y\right),(y,t)\in [0,1]\times [0,T],\hfill \\ \hfill u(1,y,t)& ={\gamma }_1{\int }_0^{1}u(x,y,t)\mathrm{d}x-2\frac{{\gamma }_1}{\pi }{\mathrm{e}}^{-t}sin\left(\pi y\right),(y,t)\in [0,1]\times [0,T],\hfill \\ \hfill u(x,0,t)& =0,u(x,1,t)=0,(x,t)\in [0,1]\times [0,T],\hfill \\ \hfill u(x,y,0)& =sin\left(\pi x\right)sin\left(\pi y\right),(x,y)\in [0,1]\times [0,1].\hfill \end{array}$$

The right-hand side function f in the differential equation and the initial and boundary conditions were prescribed to satisfy the given exact solution

$$u^{*}(x,y,t)={\mathrm{e}}^{-t}sin\left(\pi x\right)sin\left(\pi y\right)$$

of problem (1)–(5) (see, Figure 1a).

We calculated the maximum norm of the error of the numerical solution with respect to the exact solution and relative error

$$E=\underset{i=0,\cdots ,M}{max}|U_i^N-u^{*}(x_i,t_N)|,E_r=E/\left|u^{*}(x_i,t_N)\right|.$$

The results of the numerical experiment are presented in

Table 1. (Problem 1) The errors for different $h,\tau =h^2$ ($\eta =1$, ${\gamma }_0=1$, ${\gamma }_1=-1$, $T=2$).

h	$\mathit{\tau }$	E
$2^{-2}$	$6.2500·{10}^{-2}$	$4.832·{10}^{-3}$
$2^{-3}$	$1.5625·{10}^{-2}$	$1.226·{10}^{-3}$
$2^{-4}$	$3.9063·{10}^{-4}$	$3.071·{10}^{-4}$
$2^{-5}$	$9.7660·{10}^{-5}$	$7.778·{10}^{-5}$
$2^{-6}$	$2.4414·{10}^{-5}$	$1.966·{10}^{-5}$

. Note that the errors were $\mathcal{O}(\tau +h^2)$ with a sufficient accuracy for all τ and h. The ratio $\tau /h^2$ was equal to 1.

The errors E for a different ratio of step sizes τ and h are given in

Table 2. (Problem 1) The errors for different $h,\tau =h$ ($\eta =1$, ${\gamma }_0=1$, ${\gamma }_1=-1$, $T=2$).

h	$2^{-2}$	$2^{-3}$	$2^{-4}$	$2^{-5}$	$2^{-6}$
E	0.0166	0.0103	0.0056	0.0028	0.0013

. Now, $\tau /h^2=1/h$ and the error was $\mathcal{O}\left(\tau \right)$ as expected (independent of $\tau /h^2$ which varied from 4 to 64). The plots of the error distribution at mesh points for $T=2$ are presented in Figure 2a for the classical problem and in Figure 2b for the problem with nonlocal boundary conditions.

Problem 2.

In the second example, we considered the model problem (1)–(5) in the domain $D:=(0,1)\times (0,1)$. The right-hand side function f and the initial and boundary conditions were chosen so that the function

$$u^{*}(x,y,t)={\mathrm{e}}^{x+y}{\mathrm{e}}^{\alpha t+1}$$

was the exact solution of the problem (see, Figure 1b), i.e.,

$$f(x,y,t)=2\alpha {\mathrm{e}}^{x+y}{\mathrm{e}}^{\alpha t+1}\left({\mathrm{e}}^{x+y}{\mathrm{e}}^{\alpha t+1}-\eta \right)$$

$$\begin{array}{cc}\hfill u(0,y,t)& ={\gamma }_0{\int }_0^{1}u(x,y,t)\mathrm{d}x-{\gamma }_0(-{\mathrm{e}}^{\alpha t+y+1}+{\mathrm{e}}^{\alpha t+y+2})+{\mathrm{e}}^{\alpha t+y+1},(y,t)\in [0,1]\times [0,T],\hfill \\ \hfill u(1,y,t)& ={\gamma }_1{\int }_0^{1}u(x,y,t)\mathrm{d}x-{\gamma }_1(-{\mathrm{e}}^{\alpha t+y+1}+{\mathrm{e}}^{\alpha t+y+2})+{\mathrm{e}}^{\alpha t+y+2},(y,t)\in [0,1]\times [0,T],\hfill \\ \hfill u(x,0,t)& ={\mathrm{e}}^{\alpha t+x+1};,u(x,1,t)={\mathrm{e}}^{\alpha t+x+2},(x,t)\in [0,1]\times [0,T],\hfill \\ \hfill u(x,y,0)& ={\mathrm{e}}^{x+y},(x,y)\in [0,1]\times [0,1].\hfill \end{array}$$

We set $\alpha =0.1$.

In

Table 3. (Problem 2) The errors for different $h,\tau \ne h^2$ ($\eta =0.1$, ${\gamma }_0=11$, ${\gamma }_1=-10$, $T=2$).

h	$\mathit{\tau }$	E	${\mathit{E}}_{\mathit{r}}$
$0.0625$	0.005	$5.044557·{10}^{-2}$	$8.22822·{10}^{-3}$
$0.03125$	0.00125	$1.40017·{10}^{-2}$	$2.31612·{10}^{-3}$
$0.015625$	0.00031250	$3.68602·{10}^{-3}$	$6.1332·{10}^{-4}$
$0.0078125$	0.00007813	$9.4561·{10}^{-4}$	$1.5773·{10}^{-4}$
$0.00390625$	0.00001953	$2.4045·{10}^{-4}$	$4.011·{10}^{-5}$

, E tends to $\mathcal{O}(\tau +h^2)$ as τ and h decreases but is larger for large values of τ and h. It may happen because $u^{*}$ grows rapidly for large $x,y,t$. However, there might be other reasons.

Table 4. (Problem 2) The errors for different T, $\alpha =0.1$, $\eta =1$, ${\gamma }_0=11$, ${\gamma }_1=-10$, $h=0.1$, $\tau =h^2$.

T	0.5	1	2	5	10
E	0.0011	0.0018	0.0027	0.0042	0.0069

and

Table 5. (Problem 2) The errors for different $\eta$, $\alpha =0.1$, $T=2$, ${\gamma }_0=11$, ${\gamma }_1=-10$, $h=2^{-6}$, $\tau =h^2$.

$\eta$	0.1	0.5	1	2	5	10
E	$6.10·{10}^{-4}$	$5.22·{10}^{-4}$	$3.82·{10}^{-4}$	$1.63·{10}^{-4}$	$4.78·{10}^{-4}$	$8.03·{10}^{-4}$

show how the error depends on T and η. We can see that it should be proportional to the length of interval $[0,T]$ but does not depend too much on η. The plots of the error distribution for time $T=2$ are presented in Figure 3a for the classical problem and in Figure 3b for the problem with nonlocal boundary conditions.

6. Conclusions and Remarks

In the article, we investigated the stability of the ADI method for the third-order 2D linear pseudoparabolic equation with boundary integral conditions (2) and (3). The ADI method defined by Formulas (10)–(13) had not been previously studied by other authors for pseudoparabolic equations.

As is well known, one of the most important properties for any numerical method to solve differential equations with boundary and initial conditions is the stability of the method. As mentioned in the introduction, for nonlocal boundary conditions, the differential problem is not self-adjoint. Therefore, the theoretical study of the method (proof of stability and convergence) usually becomes more complicated. In other words, it is not always possible to immediately apply the well-known theoretical conclusion that convergence follows from approximation and stability. In this paper, the spectrum structure of the nonsymmetric matrix $\mathbf{S}$ was used for the theoretical study of the ADI method written in (25). In Section 4, it was proved that the stability of the differential scheme in the vector norm $\left|\right|·{\left|\right|}_{*}$ follows from the condition $\varrho \left(\mathbf{S}\right)<1$. We note that the convergence of the difference method would follow from this condition, if we proved the equivalence of the norms $\left|\right|·{\left|\right|}_{*}$ and $\left|\right|·{\left|\right|}_2$, i.e., that the inequalities $C_1\left|\right|·{\left|\right|}_{*}\le \left|\right|·{\left|\right|}_2\le C_2\left|\right|·{\left|\right|}_{*}$ were valid with constants $C_1$ and $C_2$ independent of h and $\tau$. This approach was applied in [49]. We could not claim whether that approach was valid for pseudoparabolic 2D equations with nonlocal conditions (2) and (3).

Therefore, we chose another path. As a first step, we considered a differential equation with nonlocal conditions of a new form or a new numerical method, and we chose an approach that gave us only stability (see [21,31,49,50]). Such a methodology often works when the differential equation is with constant coefficients, and the approach is related to the spectrum structure of the differential problem.

The convergence of the ADI method is in our immediate plans but requires additional research or assumptions (see Remark 3).