On the stability of explicit finite difference schemes for a pseudoparabolic equation with nonlocal conditions

Justina Jachimavičienė, Mifodijus Sapagovas, Artūras Štikonas, Olga Štikonienė Vytautas Magnus University Vileikos str. 8, LT-44404 Kaunas, Lithuania j.jachimaviciene@if.vdu.lt Institute of Mathematics and Informatics, Vilnius University Akademijos str. 4, LT-08663 Vilnius, Lithuania mifodijus.sapagovas@mii.vu.lt Faculty of Mathematics and Informatics, Vilnius University Naugarduko str. 24, LT-03225 Vilnius, Lithuania arturas.stikonas@mif.vu.lt; olga.stikoniene@mif.vu.lt


Introduction
We consider third-order linear pseudoparabolic equation with nonlocal integral conditions where η 0. The goal of this paper is to present explicit finite difference schemes (FDS) for this problem and to investigate their stability properties.
c Vilnius University, 2014 In the field of numerical analysis of differential equations, great attention has been dedicated to the development of efficient and accurate methods, mostly of implicit type.Implicit finite difference schemes for linear or nonlinear pseudoparabolic equation with Dirichlet boundary conditions were introduced in early 1970s [1,2].Stability and convergence of finite difference method were investigated in these papers.There are several situations in which the use of implicit methods causes a huge amount of computation (for example, systems of very large dimension, in presence of strong nonlinearity).Explicit methods for partial differential equations overcome these difficulties successfully.The main disadvantage of explicit methods concerns their stability.So, the benefits and stability properties of explicit methods have to be deeply analysed for safe exploiting.
Theoretical research of pseudoparabolic equations with nonlocal boundary conditions (NBC) was started due to their applications for complex problems in science and technology.One of the first results for the third-order pseudoparabolic equations with NBCs was given in monograph [3] for problems of soil dampness dynamics.A number of papers devoted to underground water flow dynamics modelling by pseudoparabolic equation with NBCs were published later, see [4,5].
Numerical methods for pseudoparabolic equations with NBCs are of permanent interest for researchers during the last decades.The Rothe time-discretization method for a nonlinear pseudoparabolic equation with integral conditions was considered in [6].A numerical method for solving a one-dimensional nonlinear pseudoparabolic equation with one integral condition in reproducing kernel space was proceeded in [7].Implicit FDS for a linear one-and two-dimensional pseudoparabolic equation with various integral conditions are analysed in [8][9][10][11].
Applications of pseudoparabolic problems with different boundary conditions including NBCs are discussed in papers [6,7] and monograph [12].Both explicit and implicit FDS for numerical solution of (1) are investigated completely enough for case of parabolic equation (η = 0).Problems with classical boundary conditions are studied in many monographs (for example, [13]).Note that, in the case η = 0, third-order derivative u txx does not allow to write two-layer explicit FDS for (1).
Three-layer FDS for the second order parabolic equation (equation (1) in the case η = 0) with classical boundary conditions are investigated in details in [13].Such schemes for problems with NBCs are studied in [14,15].Following [14], we can convert threelayer FDS to two-layer system where W j , F are vectors and matrix S is non-symmetrical due to nonlocal conditions.Energy norm of vector W j is connected with the eigenvalue problem for finite difference operator with NBCs.Nonlinear eigenvalue problem is also used for investigation of spectrum of matrix S. A number of papers devoted to the investigation of spectrum of differential or difference operators with NBCs and their applications.We refer to [16][17][18][19][20][21][22][23][24][25][26] for detailed discussions.
Applying the methodology of [14] to investigation of three-layer schemes with nonlocal conditions, we prove that stable explicit FDS can be constructed using specific approximation of third-order derivative.
The rest of this paper is organised as follows.We formulate a difference problem for approximation of the pseudoparabolic problem (1)-( 4) in Section 2. Stability regions and stability conditions are investigated in Section 3. Results of numerical experiments are presented in Section 4. Final conclusions are done in Section 5.
2 Finite difference scheme 2.1 Three-layer explicit scheme Let us write three-layer explicit FDS for differential problem (1)-(4).To this aim, we use specific approximation of the third-order derivative.As far as the authors are aware, such approximation firstly was used for solution of pseudoparabolic equation in [27,28].
We introduce grids with uniform steps We use the notation U j i := U (x i , t j ) for functions defined on the grid (or parts of this grid) ω h × ω τ .Let us define vectors U = (U 0 , U 1 , . . ., U M ) and U = (U 1 , U 2 , . . ., U M −1 ) .
If solution u of differential problem (1)-( 4) is smooth enough, then the following approximation can be written in any point (x i , t j ) ∈ ω: In order to construct explicit FDS, we approximate (u t ) j i by the forward differences, and another two derivatives by the backward differences: Let us consider three-layer FDS for problem (1)-( 4) where In order to use three-layer scheme, we need to impose an additional condition on the first layer.One way of proceeding it is to apply at the first step the two-layer implicit scheme of accuracy O(τ 2 + h 2 ) [13].
Eigenvalue problem ΛU = λU is equivalent to the eigenvalue problem for the difference operator with NBCs (see [14] for detailed explanation) Let us introduce 2(M − 1)-order vector Then three-layer FDS (13) can be written as two-layer scheme [13,14] where and 0 is zero matrix (or vector).Matrix S is nonsymmetric except case γ 0 = γ 1 = 0 (as well as matrices Λ, B and C).
In the next section, stability of FDS ( 16) will be analysed using spectrum of matrix S. So, we investigate eigenvalues of this matrix which are determined as the roots of the following equation: Equation ( 17) is characteristic equation of the nonlinear eigenvalue problem Nonlinear eigenvalue problem of such type when matrices A, B and C are symmetric is well known (see [29]).In our case, matrix B (and it is possible to investigate nonlinear eigenvalue problem (18) using another property of these matrices.Namely, matrices A, B and C have the same system of eigenvalues as matrix Λ.
Lemma 1. Eigenvalues of matrix S coincide with eigenvalues of nonlinear eigenvalue problem (18).
Lemma 2. Eigenvectors of matrices A, B and C coincide with system of eigenvectors of eigenvalue problem ( 14)-( 15) of difference operator with nonlocal conditions.
Lemma 3.All the eigenvalues of matrix Λ are real and simple.Moreover, , then all the eigenvalues are positive; , then there exists one zero eigenvalue, other eigenvalues are positive; , then there exists one negative eigenvalue, other eigenvalues are positive.
2, then all the eigenvalues of matrix Λ are nonnegative and 3 Stability of finite difference schemes

Case of real eigenvalues
Let us investigate the stability of FDS ( 6)- (10).We use matrix and vectors norms , where P is nonsingular matrix which columns are eigenvectors of matrix S and W i are coordinates of vector P −1 W. Such or similar definitions of norm for problems with NBCs are used in many works (see [15,16,19,30] and comments about norms therein).
According to definition of matrix S norm (19), stability condition for FDS ( 16) can be written as (S) 1.
Proof.Eigenvalues of matrix S can be found from equation (18).Let us denote any eigenvalue of Λ as λ and correspondent eigenvector as U.We denote eigenvalues of A, B and C as λ(A), λ(B) and λ(C), respectively.Then we have Substituting U into (18), we get quadratic equation where We rewrite this quadratic equation as with coefficients According to Hurwitz criterion, roots of ( 21) with real coefficients |µ| 1 if and only if The first inequality of Hurwitz criterion for coefficients (22) gives condition for λ The second inequality is equivalent to Condition τ λ 0 is stronger than (23).Lemma 4 implies that λ ∈ [0, β] when γ 0 +γ 1 2 and the theorem is proved.

Case of complex eigenvalues. Stability regions
In general, eigenvalues of operator Λ can be complex numbers.For example, when problem involves the integral conditions with variable weight coefficients (see ( 31)-( 32)), the eigenvalues of matrix Λ can be complex numbers [23].Such eigenvalues also can arise when the interval of integration in conditions ( 2) Only in special cases (for example, in the case of integral boundary conditions ( 15)) all the eigenvalues are real.
A polynomial satisfies the root condition if all the roots of that polynomial are in the closed unit disc of complex plane and roots of magnitude 1 are simple [32,33].
For any polynomial of the second order the following theorem [34] and lemma are valid.
Theorem 2. The roots of the second order polynomial are in the closed unit disc of complex plane and roots of magnitude 1 are simple if For quadratic equation ( 21), we have a system where λ and ϕ are unknowns.In order to obtain stability region, we find more simple necessary condition for existence of double roots instead of solving this system.If we substitute b = −2ae iϕ and c = ae 2iϕ into the first formula for coefficient b in ( 22), then we derive equation for λ From this equation we find These points are located on curve D δ in the complex plane C z (see Fig. 1).
For α = 1, equation ( 27) has two roots So, stability region is |z − 1| 1 for all ξ 0. The boundary of this region is S 1 (1), where S r (z 0 ) denotes a circle with radius r and centre z 0 .
The stability regions of a FDS for nonstationary problem (for example, parabolic, pseudoparabolic) give necessary conditions for the stability of the such FDS.If the spectrum of a corresponding difference operator belongs to the stability region, then the FDS is stable.Usually, NBCs depend on various parameters (in our case, γ 0 and γ 1 ).The spectrum of the operator Λ must lie in stability region of the investigated FDS.In general case, points of the spectrum move in the complex plane when values of the parameters in the NBCs vary.These points can be complex (problems with NBCs often are not selfadjoint).The spectrum of the operator Λ with various integral boundary conditions is investigated in [31].If we compare results about the spectrum of the operator Λ with the stability regions for the investigated FDS, then we can find parameters (and domains of them), such that FDS is stable.

Numerical experiment
We now present some numerical results demonstrating the efficiency of conditionally consistent explicit scheme.Also, we investigated the dependence of error on parameter η and time interval T .We consider a model problem ( 1)-( 4) in the case X = 1.The righthand side function f , initial and boundary conditions were prescribed to satisfy the given exact solution u * (x, t) = x(x − 1) sin t.
We consider uniform grids with different mesh sizes h and τ and analyse the convergence and accuracy of the computed solution from the present scheme.Difference problem ( 6)- (10) approximates differential problem (1)-( 4) conditionally with truncation error To demonstrate the accuracy of FDS ( 6)-( 9), we calculate the maximum norm of the error of the numerical solution as The results of the numerical test are listed for γ 0 = γ 1 = −10.
The dependence of ε versus grid step sizes τ, h and their quotient τ /h 2 is presented in Tables 1-3.
1.52588  Numerical experiment illustrates the dependence of error on the quotient τ /h 2 .We see that the numerical errors decays as the mesh size decrease.If we reduce h twice and τ eight times, then τ /h 2 reduces twice.Hence, the error of the solution also reduces, but it is difficult to evaluate how fast it reduces (see Table 1).However, this is possible for special values τ and h.If τ = h 4 , then truncation error is O(τ /h 2 + τ 2 + h 2 ) = O(h 2 ).Therefore, the error of the solution reduces 4 times as h reduces twice.These results are confirmed by simulation presented in Table 2. Similarly, for τ = h 2 , truncation error is O(1), thus the error stays roughly constant as h and τ are reduced (Table 3).
In Fig. 7 we sketch dependence of the absolute maximum errors on η for two different time intervals (T = 5 and T = 10).The error grows linearly for small η but saturates for large η and approaches some constant value.Some of the FDS for the second order parabolic equations with nonlocal integral conditions have the following property [35,36].Calculations show that the scheme is stable when T = 1, but the error becomes unbounded as T increases.Thus, FDS is not uniformly stable with respect to t. Results in Table 4(a) are important because the error of the numerical solution obtained using FDS (6)-( 9) do not grow with the interval of time.Analogous results are presented in Table 4(b) for model problem ( 1)-( 4) with exponentially growing exact solution given by u * (x, t) = x 3 exp t.
5 Conclusions and remarks FDS ( 6)-(10) based on ( 5) is only conditionally consistent.It means that it is not sufficient to reduce step sizes h and τ to improve the accuracy of the solution.It is necessary to reduce τ /h 2 as well.This is fairly strong restriction.However, our attempts to approximate derivative u txx better resulted in unstable FDS.
The FDS with equation instead of equation ( 6) is unstable as well.In this case, the truncation error is O(τ + h 2 ).
In order to get a conditionally stable scheme, we can leave the derivative u txx approximated by (5) while two other derivatives u t and u xx are approximated in a different way than equation (28).
Numerical experiment allows us to conclude that stable explicit schemes with conditional approximation are useful in practice.Even though step size τ can be fairly small, it is not a serious problem for modern computers.Explicit methods are very suitable for parallel computing unlike the implicit schemes.
Let us note another fact about the stability of FDS for pseudoparabolic equations.The analysis presented in this work is independent on boundary conditions.It is important that the boundary conditions must be such that eigenvalues of matrix Λ would be nonnegative.For example, we can take more general boundary conditions with weights u(0, t) = X 0 α(x)u(x, t) dx + v l (t), 0 t T, u(X, t) = X 0 β(x)u(x, t) dx + v r (t), 0 t T, instead of nonlocal conditions (2), (3).Restrictions for coefficients α(x) and β(x) when Λ will have only positive eigenvalues are investigated in [22,23,30,36].In the case of complex eigenvalues, we can investigate the stability using the method described in Subsection 3.2.