In numerical solving time-dependent boundary value problems, the emphasis is on using computational algorithms of higher accuracy (see, e.g., HundsdorferVerwer2003 ; Gustafsson2008 ). Along with increasing the accuracy of discretization in space, increasing the accuracy of discretization in time is also considered taking into account numerical methods developed for ordinary differential equations Ascher ; LeVeque
. In view of features of unsteady problems for partial differential equations, we, first of all, should be guided by methods of numerical solving Cauchy problems for stiff systems of ordinary differential equationsHairerWanner2010 ; Butcher2008 .
As a rule, discretization in time for numerical solving boundary value problems for parabolic equations is constructed on the basis of two- or three-level schemes. To study stability of such difference schemes, the theory of stability and well-posedness for operator-difference schemes Samarskii1989 ; SamarskiiMatusVabishchevich is employed. In particular, stability conditions for standard schemes with weights (-method) are obtained for linear problems in different Hilbert spaces. For two-level schemes, where the solution at two consecutive time levels is involved, polynomial approximations are explicitly or implicitly employed for operators of difference schemes. The Runge-Kutta methods (see Butcher2008 ; DekkerVerwer1984 ) are well-known examples of such schemes widely used in modern computational practice. The main feature of multilevel schemes (multistep methods) results in the approximation of time derivatives with higher accuracy on a multipoint stencil. Multistep methods based on numerical backward differentiation formulas Gear1971 should be mentioned as a typical example.
In some cases, it is possible to construct exact difference schemes. The main possibilities of their construction for solving Cauchy problems for ordinary differential equations are reflected in the works mickens1994nonstandard ; mickens2002nonstandard ; anguelov2001contributions . The idea of increasing accuracy by modifying approximations of time derivatives can be realized not only for linear equations, but also for nonlinear ones and for systems of equations. Separately, we highlight the works mickens1999nonstandard ; anguelov2005non ; patidar2005use ; hernandez2013nonstandard , where such a technology of constructing discretization in time is applied to solving parabolic equations. Using perturbations with a low parameter for standard difference approximations, we can provide a new quality of the solution, which is associated, for example, with the fulfillment of some conservation law.
For a Cauchy problem for a homogeneous linear parabolic equation, the solution can be represented as a superposition of modes, which are associated with the corresponding eigenfunctions of the elliptic operator. Dominant modes are the most important, since higher modes decay very quickly. Such a qualitative behavior of the solution must hold for a discrete problem as well. For long times, the fundamental mode of the solution is highlighted, when a regular regime is considered (see, e.g.,luikov2012analytical ). In the theory of difference schemes samarskii1996computational , there is recognized a class of asymptotically stable difference schemes that provide the correct behavior of the approximate solution for large times. SM (Spectral Mimetic) properties of operator-difference schemes for numerical solving the Cauchy problem for evolutionary equations of first order are associated with the time-evolution of individual modes of the solution. For long times, the fundamental mode of the solution is highlighted for a regular regime. With this in mind, SM-stable difference schemes are constructed Vabischevich2010b .
Here we construct difference schemes for parabolic equations, which are exact for the most important component of the solution, namely, the fundamental mode is calculated exactly. Two-level fundamental mode exact schemes (FMES) are constructed using Padé approximations.
The paper is organized as follows. A boundary value problem for a linear parabolic equation with a self-adjoint elliptic operator of second order is considered in Section 2. Using finite element approximations in space, we formulate a Cauchy problem for the corresponding differential-operator equation. In Section 3, we consider standard two-level schemes with weights and formulate stability conditions. Schemes that are exact for the solution associated with the fundamental eigenfunction of the corresponding elliptic operator are constructed in Section 4. Section 5 presents the results of numerical experiments on studying FMES accuracy for a model initial-boundary value parabolic problem.
2 Problem statement
Let be a bounded domain () with a piecewise smooth boundary . Define an elliptic operator such that
The operator is defined on the set of functions that satisfy on the boundary the following conditions:
Coefficients , and are smooth functions in and
We consider the Cauchy problem
where, for example, with the notation .
Let be the scalar product and norm in , respectively:
Multiplying equation (3) by and integrating over the domain , we arrive at the equality
Here is the following bilinear form:
with a constant . In view of (4), we put
), we have the a priori estimate
To show this, we put in (5) and get
Taking into account the lower bound for the bilinear form , we get
For numerical solving the initial-boundary value problem (3), (4), discretization in space is constructed using the finite element method brenner ; Thomee2006 . The weak formulation (5), (6) is employed. Define the subspace of finite elements and the discrete elliptic operator as
The operator acts on the finite dimensional space and
where is the identity operator.
where with denoting -projection onto . Similarly to (7), we establish the estimate
we have real eigenvalues
The eigenfunctions form a basis in . Therefore, for , we have
With this representation, we get the a priori estimate (1) with .
Various effects arising from modes can be taken into account in constructing difference schemes. Assume that the fundamental mode is separated from other modes, i.e., . Then, for large times, a regular regime is established for the solution:
If such a property is preserved after discretization in time, then we speak of asymptotic stability of difference schemes samarskii1996computational . In the representation (12), higher modes decay faster with time. Such a monotone behavior of different modes is inherited in SM-stable difference schemes Vabischevich2010b .
It seems natural to construct discretizations in time that are more accurate in reproducing the time dependence of dominant modes. In particular, we distinguish nonstandard difference schemes, which are exact for the fundamental mode of the solution, namely, fundamental mode exact schemes.
3 Nonstandard schemes with weights
The standard approach to the solution of the problem (9), (10) is associated with using two-level schemes. Let be a step of a uniform grid in time such that , . For a constant weight parameter , we approximate equation (9) by the following two-level scheme:
with the initial conditions
For any , the equation
It is possible to select a weight value in such a way that the fundamental mode is calculated exactly. The value such that is shown in Fig. 1. We restrict ourselves to the values . In the most important case, we have . For computational practice, the fully implicit scheme () is more interesting, but it cannot describe accurately the time-variation of the fundamental mode.
In our previous work afanas2013unconditionally , special unconditionally stable schemes were constructed for convection-diffusion problems using a transition to new unknowns. The main feature of these schemes is connected with the negativity of the constant in the lower bound for the operator (the estimate (8) in the present work). Here we use this idea for constructing fundamental mode exact schemes.
In view of (17), suppose that . Thus, we arrive at the scheme
with the initial condition
A comparison of the nonstandard scheme (20), (23), (24) with the standard scheme with weights (13), (14) indicates that we modified not only discretization in time, but also the problem operator itself:
Using the representation (15), we obtain
Thus, the condition (16) is fulfilled and this scheme is of FMES-type.
The nonstandard scheme (20), (23), (24) is unconditionally stable under the standard restriction on the weight . The proof is based on the well-known (see, e.g., Samarskii1989 ; SamarskiiMatusVabishchevich ) condition for the stability of the scheme with weights (21), (22). Multiply equation (18) scalarly by
In view of , we get
Thus, for , we have
4 Schemes based on Padé approximations
Two-level finite difference schemes of higher approximation order for time-dependent linear problems can be conveniently constructed on the basis of Padé approximations for the corresponding operator (matrix) exponential function Higham2008 . In the case of linear systems of ODEs, such approximations correspond to various variants of the Runge-Kutta method HairerWanner2010 ; Butcher2008 ; DekkerVerwer1984 .
For the corresponding difference scheme, we have
It is necessary to choose one or another function , which approximates the function for, where . To obtain an exact scheme for the fundamental mode, we must (see (16)) fulfill .
The Padé approximation of the function is
where and are polynomials of degrees and , respectively. These polynomials have BakerGraves-Morris1996 the form
At the point the approximation is exact, i.e. . We search the exact approximation at the point . Taking this into account, we put
In fact, this means that the transformation (25) is realized. Thus, we can formulate the following statement.
Let us highlight the most interesting variants of the constructed schemes of FMES-type. The scheme (27), (29) with corresponds to the fully implicit scheme ( in(23), (24)). The scheme (27), (29) with is associated with the Crank-Nicolson scheme ( in (23), (24)).
Only the schemes with , i.e.,
where the function
is a truncated Taylor series for .
In the class of schemes (30), in addition to the fully implicit scheme of the first-order approximation in time (), we highlight the scheme with . In this case, the approximate solution is determined from
5 Numerical experiments
The model problem (1)–(4) is considered in the unit square (). We assume that , and the coefficient in equation (1) is discontinuous. Namely, in the left lower quarter of the computational domain , it is equal to 10, in other parts of the region it equals 1. Neumann-type or Robin-type boundary conditions are specified on different parts of the boundary. The values of the coefficient and function in the boundary condition (2) are shown in Fig. 2. Thus, we consider the problem (3), (4), where
A uniform grid in space is used. Piecewise-linear finite elements on triangles are employed.
The constructed nonstandard difference schemes are based on using the known fundamental eigenvalue . The computation of the first eigenvalues and the corresponding eigenfunctions is the standard problem of computational mathematics Saadbook . We can focus on applying well-designed algorithms and relevant open-source software. To solve the spectral problems, we recommend the SLEPc package (Scalable Library for Eigenvalue Problem Computations, http://slepc.upv.es/). In this library, in particular, there are implemented the Krylov-Schur algorithm and a variant of the Arnoldi method proposed by stewart2002krylov .
We apply, for example, the standard inverse iteration method to find the fundamental eigenvalue of the spectral problem:
In the simplest case, we can put as the initial value. At the -st iteration, we have
For the -st approximate values of the first eigenvalue and the corresponding eigenfunction, we get
For the operator , we have the same eigenfunctions, and . The fast convergence of the inverse iteration method for the simplest variant (32) is presented in Table 1. The calculations are performed on the sequence of grids derived via dividing the size of triangles by half. A small number of iterations is enough to find the first eigenvalue and the corresponding eigenfunction. The first eigenfunction on the grid with is shown in Fig. 3. The solution value on contour levels is a multiple of .
The time-dependent problem is considered until . The initial condition (4) is taken in the form . For the basic variant, we put . We compare the standard fully implicit scheme (( in (13)) and the nonstandard scheme (23) with the same value . The accuracy of calculating the amplitude of the fundamental mode is estimated by the following value:
Numerical results for the amplitude of the fundamental mode obtained on different grids in time are given in Fig. 4.
It is possible to evaluate the general error of the approximate solution without separating the fundamental mode. As a reference solution (the benchmark solution used for a comparison) we employ the solution that is obtained using the fully implicit scheme (13), (14) on the finest grid in time (the number of time steps is ). The accuracy of the approximate solution is estimated at each time level:
Numerical results for the fully implicit scheme (13), (14) are given in Fig. 5. Similar data for the nonstandard scheme (23), (24) for are shown in Fig. 6. We observe a much higher accuracy of the approximate solution obtained using the FMES-type scheme. The influence of the constant is illustrated in Figs. 7 and 8. Obviously, the nonstandard scheme demonstrates higher accuracy in comparison with the standard scheme.
The publication was financially supported by the Ministry of Education and Science of the Russian Federation (the Agreement # 02.a03.21.0008).
- (1) W. H. Hundsdorfer, J. G. Verwer, Numerical Solution of Time-dependent Advection-diffusion-reaction Equations, Springer Verlag, 2003.
- (2) B. Gustafsson, High Order Difference Methods for Time Dependent PDE, Springer Verlag, 2008.
- (3) U. M. Ascher, Numerical Methods for Evolutionary Differential Equations, Society for Industrial Mathematics, 2008.
- (4) R. J. LeVeque, Finite Difference Methods for Ordinary and Partial Differential Equations: Steady-state and Time-dependent Problems, Society for Industrial Mathematics, 2007.
- (5) E. Hairer, G. Wanner, Solving Ordinary Differential Equations II: Stiff and Differential-algebraic Problems, Springer Verlag, 2010.
- (6) J. C. Butcher, Numerical Methods for Ordinary Differential Equations, Wiley, 2008.
- (7) A. A. Samarskii, The Theory of Difference Schemes, Marcel Dekker, New York, 2001.
- (8) A. A. Samarskii, P. P. Matus, P. N. Vabishchevich, Difference Schemes with Operator Factors, Kluwer Academic Publishers, 2002.
- (9) K. Dekker, J. G. Verwer, Stability of Runge-Kutta Methods for Stiff Nonlinear Differential Equations, North-Holland, Amsterdam, 1984.
- (10) C. W. Gear, Numerical Initial Value Problems in Ordinary Differential Equations, Prentice Hall, NJ, 1971.
- (11) R. E. Mickens, Nonstandard Finite Difference Models of Differential Equations, world scientific, 1994.
- (12) R. E. Mickens, Nonstandard finite difference schemes for differential equations, The Journal of Difference Equations and Applications 8 (9) (2002) 823–847.
- (13) R. Anguelov, J. M.-S. Lubuma, Contributions to the mathematics of the nonstandard finite difference method and applications, Numerical Methods for Partial Differential Equations 17 (5) (2001) 518–543.
- (14) R. E. Mickens, Nonstandard finite difference schemes for reaction-diffusion equations, Numerical Methods for Partial Differential Equations 15 (2) (1999) 201–214.
- (15) R. Anguelov, P. Kama, J.-S. Lubuma, On non-standard finite difference models of reaction–diffusion equations, Journal of computational and Applied Mathematics 175 (1) (2005) 11–29.
- (16) K. C. Patidar, On the use of nonstandard finite difference methods, Journal of Difference Equations and Applications 11 (8) (2005) 735–758.
- (17) E. Hernandez-Martinez, H. Puebla, F. Valdes-Parada, J. Alvarez-Ramirez, Nonstandard finite difference schemes based on Green’s function formulations for reaction–diffusion–convection systems, Chemical Engineering Science 94 (2013) 245–255.
- (18) A. Luikov, Analytical Heat Diffusion Theory, Academic Press, 1968.
- (19) A. A. Samarskii, P. N. Vabishchevich, Computational Heat Transfer, Wiley, 1996.
- (20) P. N. Vabishchevich, Two-level finite difference scheme of improved accuracy order for time-dependent problems of mathematical physics, Computational Mathematics and Mathematical Physics 50 (1) (2010) 112–123.
- (21) S. C. Brenner, R. Scott, The Mathematical Theory of Finite Element Methods, Springer, 2008.
- (22) V. Thomée, Galerkin Finite Element Methods for Parabolic Problems, Springer Verlag, 2006.
- (23) N. M. Afanas’eva, P. N. Vabishchevich, M. V. Vasil’eva, Unconditionally stable schemes for convection-diffusion problems, Russian Mathematics 57 (3) (2013) 1–11.
- (24) N. J. Higham, Functions of Matrices: Theory and Computation, SIAM, Philadelphia, 2008.
- (25) G. A. Baker, P. R. Graves-Morris, Padé Approximants, Cambridge University Press, 1996.
- (26) Y. Saad, Numerical Methods for Large Eigenvalue Problems, SIAM, 2011.
- (27) G. W. Stewart, A Krylov–Schur algorithm for large eigenproblems, SIAM Journal on Matrix Analysis and Applications 23 (3) (2001) 601–614.