On Construction and Analysis of Finite Difference Schemes for Pseudoparabolic Problems with Nonlocal Boundary Conditions

AbstractIn this paper the one- and two-dimensional pseudoparabolic equations with nonlocal boundary conditions are approximated by the Euler finite difference scheme. In the case of classical boundary conditions the stability of all schemes is investigated by the spectral method. Stability regions of finite difference schemes approximating pseudoparabolic problem are compared with the stability regions of the classical discrete parabolic problem. These results are generalized for problems with nonlocal boundary conditions if a matrix of the finite difference scheme can be diagonalized. For the two-dimensional problem an efficient algorithm is constructed, which is based on the combination of the FFT method and the factorization algorithm. General stability results, known for the three level finite difference schemes, are applied to investigate the stability of some explicit approximations of the two-dimensional pseudoparabolic problem with classical boundary conditions. A connection between the energy met...


Introduction
The correctness of main types of boundary value problems for parabolic equations is well investigated in many papers and textbooks. It is shown that boundary conditions can introduce essential changes in construction and theoretical analysis of stable numerical approximations of differential problems with nonlocal boundary conditions. In various real-world applications non-classical boundary conditions are formulated and nonlocal boundary value problems make a basis for important modelling and simulation projects. Thus the analysis of mathematical models and numerical methods for solving such problems is considered in many papers, e.g. parabolic problems in heat conduction and thermodynamics are investigated in [1,3,8], pseudoparabolic problems in underground water flow are studied in [2,9,28].
Let us consider domain D = (0, 1). We formulate an initial boundary value pseudoparabolic problem: where f , ϕ, k, g are known functions, such that We also consider nonlocal boundary conditions: u(x, t) dx, t > 0, (1.5) where γ 1 , γ 2 are given constants. The stability analysis of numerical approximations of parabolic and pseudoparabolic problems with nonlocal boundary conditions is done by various methods. Here we review the three main classes of such methods.
For parabolic problems with nonlocal boundary conditions, the maximum principle is the most popular technique for the stability analysis, see [7,10]. The method of bounding functions is used in [8] to get the stability estimate in the maximum norm. We note, that for stationary problems this method enables us to get sufficient (and in many cases even necessary and sufficient) stability conditions (see, [6] and references contained therein).
The energy method is used to investigate the stability of the numerical approximations in the case of general parabolic and pseudoparabolic equations with non-constant coefficients, see [2,10,12,14] and references contained therein. The drawback of this method is that in order to apply energy estimates we should answer the important question how to find the right norm or functional suited for the given problem. In this paper we apply general techniques which are known for the two-level and three-level operator based finite difference or finite volume approximations of partial difference equations [23]. Then some simple sufficient stability conditions are known for the given templates of operator equations.
A general technique to prove necessary and sufficient stability conditions for non-stationary numerical approximations is the spectrum method. Our aim is to show a connection of the energy based techniques and spectral stability analysis methods. It is important to note, that for non-normal matrices (which are typical in the case of nonlocal boundary conditions) a modified spectral method technique gives a similar eigenvalue criterion [16,17]. This similarity makes a bridge for application of stability results known for pseudoparabolic problems with classical boundary conditions to problems with nonlocal boundary conditions. The main features of the spectral technique for non-normal operators can be explained for a linear system where λ j are the eigenvalues of A. Let us introduce vectors Z(t) = U −1 w(t), then: If we know that max j Re λ j ≤ ω, then the stability estimate follows: In general, if the matrix A is not normal, an a priori estimate of cond (U ) may be difficult to obtain in a suitable norm. This technique is applied for the stability analysis of parabolic and pseudoparabolic problems with nonlocal boundary conditions in [15,18,19,20]. We note, that in this paper we will not investigate the spectrum of obtained discrete schemes with nonlocal boundary conditions. The stability analysis of the constructed finite difference schemes is based on the assumption that a matrix of the discrete operator can be diagonalized and eigenvectors make a complete basis system. Then the obtained results on the stability regions of discrete approximations of pseudoparabolic problems can be used also in the case of nonlocal boundary conditions. Here we also note that the structure of the spectrum of some differential operators with nonlocal boundary conditions is investigated in [25,26].
Starting from pioneer works [11,13], many numerical algorithms have been developed for equations of pseudoparabolic type with various classical boundary conditions. The explicit-implicit Euler and Crank-Nicolson approximations are constructed in [13]. The ρ-stability and convergence analysis is done in the discrete L 2 norm by using the spectral analysis technique. Some comparison of approximations for parabolic and pseudoparabolic equations is presented, and it is noted that even the backward time pseudoparabolic problems can be solved by the proposed discrete algorithms.
Our first goal is to define the stability regions of two-level finite difference schemes. A class of finite difference schemes with averaging in time is considered, including high-order approximations, stability regions of these schemes are constructed and investigated. Then by using a general template defined above, the obtained results on the stability regions of the discrete pseudoparabolic operators are applied for problems with nonlocal boundary conditions. The essential assumption is used that matrices of the constructed finite difference operators can be diagonalized and eigenvectors define a complete basis set.
The second goal of this paper is to investigate the stability of three-level explicit approximations of multidimensional pseudoparabolic problems. There are many papers where implicit finite difference and finite element approximations are used, see [11,22]. In order to find a discrete solution at each time level a large system of linear equations is solved by iterative methods. There we consider two approaches to construct efficient numerical algorithms to solve multidimensional pseudoparabolic problems with various boundary conditions. The first one is to use explicit finite difference schemes. Here we analyze the stability of explicit algorithms, that are constructed as multidimensional generalizations of one-dimensional algorithms proposed in [18]. Our goal is to apply general results of the stability of three-level finite difference schemes. First we consider the case of classical boundary conditions, when discrete operators are symmetrical and positive. It is shown that the technique of operator estimates enables users to derive the stability conditions for all proposed approximations. Next we show a simple connection of the operator based estimates with standard spectral estimates based on the root condition of characteristic equations (again for symmetrical and positive operators). This result enables us to investigate the stability of explicit schemes for multidimensional pseudoparabolic problems with nonlocal conditions. These stability results are valid if the assumption on diagonalization of a matrix of the discrete operator holds and all eigenvalues are real.
The second possibility is to use splitting type approximations. A review of this approach is given in [27]. Our aim is to investigate the convergence of one Locally One-Dimensional (LOD) scheme, which is proposed in [19]. The stability of this scheme is investigated in [19] by using the same assumption that a matrix of the discrete operator can be diagonalized and the eigenvectors define a complete basis set. Thus the main task is to investigate the approximation error the discrete scheme in appropriate norms.
The rest of this paper is organized as follows. The structure of the spectrum of the pseudoparabolic equation is investigated in Section 2. In Section 3, the implicit Euler finite difference scheme is constructed for approximation of the pseudoparabolic problem (1.1)-(1.3) and the accuracy of the approximation is investigated. The stability regions of the finite difference schemes are derived and investigated in Section 4. In Section 5, the two-dimensional pseudoparabolic problem with nonlocal boundary conditions is solved. The approximation and convergence of the split step and the full approximation implicit finite difference schemes is investigated. Explicit three-level multidimensional finite difference schemes are constructed in Section 6. Some final conclusions are done in Section 7.

The Spectral Stability Analysis
Let us write the differential equation (1.1) in the operator form where It is easy to check, that in the case of boundary conditions (1.3) these operators are symmetric in the Hilbert space where (u, v) denotes the inner product of functions u and w: With this notation the eigenvalue problem can be stated as: find a number λ and a function ϕ ∈ H 1 0 (D), ϕ = 0 such that Let λ j be the jth eigenvalue and ϕ j the corresponding eigenfunction. We assume that eigenfunctions {ϕ j } ∞ j=1 form an orthonormal basis for H 1 0 (D). This statement is true if operators A and B are symmetric and positive definite.
We now seek the solution to problem (1. Let us assume that f ≡ 0. Inserting the sum (2.3) into the differential equation (1.1) and using (2.2) we obtain the ODE for functions c j : We obtain the standard conclusion, that the given differential problem is stable if Re λ j ≥ 0, but for the pseudoparabolic problem the eigenvalues are defined by the generalized eigenvalue problem (2.2) (see also [13]). Next we consider the important example, when Let us assume that eigenvalues and eigenvectors of the operator A are defined as Aϕ j = µ j ϕ j , j = 1, 2, . . . .
Then we may easily determine the eigenvalues of the generalized problem (2.2) explicitly For the definition of the stability region, we consider a general case of complex eigenvalues. Let us denote µ = µ R + iµ I . After simple computations we find that the stability region of the pseudoparabolic problem is defined by the inequality By assuming that η > 0, we can write it in the following form Note, that the obtained stability region is much larger than the stability region µ R ≥ 0 of the classical parabolic equation.

Finite Difference Schemes
For the construction and analysis of finite difference schemes we shall use results obtained in [8]. The domainD is covered by the discrete uniform grid where τ is the time step. Although the constant time step is taken here, the following studies can be easily extended to the case when τ varies. We consider numerical approximations U n j to the exact solution values u n j = u(x j , t n ) at the grid points (x j , t n ) ∈D h × ω τ . For functions defined on the grid we introduce the forward and backward difference quotients with respect to x and similarly the backward difference quotient and the averaging operator with respect to t We approximate the differential problem (1.1)-(1.3) by the finite difference scheme Let us define discrete operators  2) approximates the differential problem and the following estimate is valid for the residual error where the residual is defined by Proof. It is sufficient to analyze the approximation error of the term B h ∂tu n , the remaining analysis can be done in a similar way. Let us denotet := t n−1/2 . Using the Taylor formula with the remainder term in the integral form we obtain Then we obtain the estimate Let us consider the important example [20] k(x) = 1, Then we construct the high-order approximation finite-difference scheme where the notation LU j := ∂x∂ x U j is used. This scheme is defined on the same compact stencil of the grid ω h as for the finite difference scheme (3.1) and the approximation accuracy of (3.5) is equal to O(τ 2 + h 4 ).

The Stability Region
In this section we find the stability regions of the forward and backward Euler methods. Here we restrict ourselves to the case B = ηA.
The forward Euler scheme. We take θ = 0 in the equation (3.1). The scalar stability test equation is defined by where µ = µ R + iµ I are complex eigenvalues of the operator A. Then we get the recursion The stability region of the finite difference scheme is After simple computations we obtain that the stability condition for the forward Euler scheme is defined by If τ < 2η, then the stability region of the forward Euler finite difference scheme is The backward Euler scheme. We take θ = 1 in the equation (3.1). The scalar stability test equation is defined by Then we get the recursion c n = qc n−1 , q = 1 + ηµ 1 + ηµ + τ µ .
Real eigenvalues. Now we consider the case when operators A and B are symmetrical and positive definite Then the stability analysis can be based on general results of the stability of two-level finite difference schemes [23]. Let us consider the problem where I + τ R > 0 and the operator R is called a regularizator. If A = A * > 0, then a sufficient stability condition is given by For the pseudoparabolic discrete problem (3.1) the regularizator is equal to Let us assume that B ≥ cA, then the finite difference scheme (3.1) is stable if the inequality is valid. It follows from (4.1) that the symmetrical scheme with θ = 1 2 and the high-order approximation scheme θ = 1 2 − h 2 12τ are unconditionally stable.
Nonlocal boundary conditions. In this section, the stability estimates are derived for the solution of problem (3.1)-(3.2), when the discrete nonlocal boundary conditions are formulated instead of boundary conditions (3.3). We approximate the integral conditions (1.5) by the discrete conditions U n 0 = γ 1 S J U n , U n J = γ 2 S J U n , n > 0, (4.2) where S J U n approximates the integral by the trapezoidal rule The obtained finite difference scheme (3.1), (3.2), (4.2) defines the nonsymmetrical operators A, B. Let us consider the case of operators A = −L, B = −ηL. Then it is proved in [24] that for γ 1 + γ 2 < 2 all eigenvalues of the operator A are real and positive. Then the stability of the finite difference scheme (3.1), (3.2) can be investigated by using the matrix diagonalization method (1.6). The simple spectral stability analysis proves that the stability condition (4.1) is valid for the problem with nonlocal boundary conditions, if we use the norm A = max j |µ j |. This matrix norm is compatible with the vector norm X = U −1 X ∞ , where U is the matrix of linearly independent eigenvectors of the matrix A.
The main goal of this section is to investigate efficient numerical approximations of the multidimensional pseudoparabolic problem. Parallel numerical algorithms for the three-dimensional parabolic problem with nonlocal conditions are investigated in [4].
The LOD finite difference scheme. The following splitting into locally one-dimensional (LOD) differential problems is proposed in [19]. This type of splitting is based on ideas used for solving the multidimensional parabolic problems [23]. It is proved in [19] that the computational algorithm which is obtained after approximation of the one-dimensional sub-problems by the finite difference schemes is stable. This result is obtained by using the diagonalization of the matrix of two-dimensional discrete elliptic operator with discrete nonlocal boundary conditions. We note here that such a stability result is not sufficient to prove the convergence of the discrete solution of the LOD scheme to the solution of the differential problem (5.1)-(5.5). We will show that this LOD scheme is not approximating the two-dimensional problem (5.1)-(5.5). In order to prove this result, let us consider the following scalar test problem: where a 1 , a 2 > 0 are given positive constants. The exact analytical solution of equation (5.8) is Next, let us consider the split step approximation of (5.8) The solution of this discrete problem is equal to We see that the error Z n = C n − c n does not depend on the split parameter τ and therefore the convergence can not be proved for the solution of the LOD scheme (5.9).
The backward Euler finite difference scheme. In this section we construct the efficient algorithm to solve the two-dimensional pseudoparabolic problem. It is based on spectral method and FFT. The domainD is covered by the discrete uniform grid D h = (x j , y k ): x j = jh, y k = kh, j, k = 0, . . . , J , Jh = 1.
We take into account the structure of the boundary conditions (5.13) and introduce the Fourier sums ϕ n α,k sin πky j .
Substituting these sums into finite difference scheme (5.10)-(5.14) we obtain a sequence of independent discrete problems for spectral coefficients U n k,i , k = 1, . . . , J − 1: U n k,0 = γ 1 S h U n k + ϕ n 1,k , U n k,J = γ 2 S h U n k + ϕ n 2,k . Let us assume that the parameters γ 1 and γ 2 are such (or the time step τ is sufficiently small) that the obtained systems of linear equations are solvable. Then all these problems can be solved independently in parallel by using the modified factorization algorithm tailored for solving one dimensional discrete problems with the nonlocal boundary conditions [8]. The complexity of this part of the algorithm is O(J 2 ) floating point operations.
The Fourier sums are computed by using the FFT algorithm. If J = 2 k the complexity of this part of the scheme is O(J 2 log J) operations. Thus the full complexity of the implementation of one time step of the finite difference scheme (5.10)-(5.10) is O(J 2 log J).
The presented implementation algorithm enables us also to prove the stability of the given finite difference scheme, if the assumption that matrices of one-dimensional operators can be diagonalized and the eigenvectors define a complete set of vectors. Note that results given in Section 5 on the stability of 1D problems can be used here.

Explicit Finite Difference Schemes
In this section we consider the stability of some explicit finite difference schemes, which were proposed and analyzed in [18]. We will generalize these schemes for two-dimensional pseudoparabolic problem and note that three-dimensional problems can be solved in a similar way. Explicit schemes are very important since parallelization of these algorithms can be done very efficiently.
For the most part of our analysis we consider classical boundary conditions. Our goal is to investigate the efficiency of a general stability analysis technique which is developed for three-level operator equations [23]. Also, we are interested to establish the connection between this energy method and the Hurwitz stability criterion, since the last one can be used for problems with nonlocal boundary conditions. Let us introduce the difference operators The well-known stability result is valid for the three-level finite difference schemes written in a canonical form [23]: Theorem 1. Let operators R and A be symmetric and positively definite A = A * > 0, R = R * > 0 and they do not depend on t. Then the following conditions B(t) ≥ 0, R > 1 4 A are sufficient for the stability of (6.1) with respect to the initial conditions. This stability criteria was applied for the analysis of hyperbolic heat conduction equations in [5,21].
Finite difference scheme 1. We solve the two-dimensional equation (5.1). The most simple explicit three-level finite difference scheme is obtained if both diffusion operators are approximated on the nth time level: The initial condition for n = 1 time level can be computed by using the implicit scheme (3.1). This selection has no influence on the stability of the explicit scheme (6.2). It is easy to check that the approximation error of this explicit finite difference scheme is O(τ + h 2 ). The finite difference scheme (6.2) can be written in the canonical form by taking Then we straightforwardly obtain, that the stability condition R > 1 4 A is not satisfied, therefore the scheme (6.2) is unstable.
Finite difference scheme 2. The second scheme is constructed as a modification of the well-known Dufort-Frankel three-level finite difference scheme [23]: The approximation error of this scheme is O(τ + h 2 + τ h 2 ). Finite difference scheme (6.3) can be written in the canonical form by taking From A < 8 h 2 I it follows that the Dufort-Frankel scheme is unconditionally stable for the two-dimensional parabolic problem when the parameter η = 0. But since (AW , W ) > 4 h 2 (IW , W ) for some eigenvectors W , the explicit finite difference scheme (6.3) is still unstable.
Finite difference scheme 3. The third explicit finite difference scheme is constructed by perturbing the finite difference scheme (6.2): The approximation error of this scheme is O(τ + h 2 + τ h 2 ). Finite difference scheme (6.4) can be written in the canonical form by taking It follows from the definition of the operators A and R that the inequality R > 1 4 A is satisfied if the restriction τ ≤ h 2 is valid. The obtained restriction on discrete step τ and h improves the results of [18], where the conditional stability estimate is proved for τ < h 2 4 (in the case of one-dimensional problem).
Finite difference scheme 4. Let us consider the following explicit finite difference scheme The accuracy of this approximation again is of order O(τ + h 2 + τ h 2 ). In [18], one-dimensional version of scheme (6.5) was used in numerical experiments and no stability results are presented. It can be written in the canonical form by taking Applying the stability criteria of Theorem 1 we obtain, that the explicit threelevel scheme (6.5) is unconditionally stable.
Nonlocal boundary conditions. As was noted above, in the case of the nonlocal boundary conditions the Neumann stability analysis is changed to the analysis of the spectrum of non-symmetrical operators. First, we restrict ourselves to the cases when the diagonalization of the matrix can be done and all eigenvalues are real. In the case of the three-level finite difference schemes the stability of the difference scheme depends on the solutions of the quadratic characteristic equation b q 2 − 1 2τ + r q 2 − 2q + 1 + aq = 0, (6.6) where a, b, r are given constants. The well known Hurwitz stability criterion defines that the solution of the characteristic equation (6.6) satisfies the estimate |q| < 1 if the following estimates are valid: After simple computations we obtain the equivalent conditions b > 0, a > 0, a 4 < r, (6.7) which are the scalar variants of the energy estimates from Theorem 1. Thus we can apply conditions (6.7) for the stability analysis of explicit finite difference schemes in the case of nonlocal boundary conditions. Next we examine shortly a case when eigenvalues of the matrix are complex. There is no direct generalizations of the Hurwitz stability criterion for complex quadratic polynomials to give us simple constructive conditions of the stability. One possible way how to investigate the stability of explicit finite difference schemes is to construct a contour of the stability region by using the techniques developed for general multistep numerical algorithms [17]. But for pseudoparabolic problems the contour depends on parameters τ , h and η. Another possibility is to derive sufficient stability conditions from some wellknown bounds on the roots of a characteristic polynomial, e.g. the Rouche theorem can be used.

Conclusions
In this paper the one-and two-dimensional pseudoparabolic equations with nonlocal boundary conditions are approximated by the Euler finite difference scheme. It is proved that the stability regions of the schemes for pseudoparabolic problem are much larger than for the classical parabolic problem. For the two-dimensional problem the efficient algorithm is constructed, which is based on the combination of the FFT method and the factorization algorithm. It is proved that the general stability results of the three level finite difference schemes can be applied efficiently to investigate the stability of explicit approximations of the pseudoparabolic problem. The unconditionally stable explicit schemes are proposed. a qualitative connection between the stability conditions derived by the energy method and the spectrum Hurwitz stability criterion is shown, which enables to transfer directly the obtained stability results also for pseudoparabolic problems with nonlocal boundary conditions if the matrix of the discrete operator can be diagonalized and all eigenvalues are real numbers.