1 Introduction
In numerical solving timedependent 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 equations
HairerWanner2010 ; Butcher2008 .As a rule, discretization in time for numerical solving boundary value problems for parabolic equations is constructed on the basis of two or threelevel schemes. To study stability of such difference schemes, the theory of stability and wellposedness for operatordifference 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 twolevel schemes, where the solution at two consecutive time levels is involved, polynomial approximations are explicitly or implicitly employed for operators of difference schemes. The RungeKutta methods (see Butcher2008 ; DekkerVerwer1984 ) are wellknown 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 operatordifference schemes for numerical solving the Cauchy problem for evolutionary equations of first order are associated with the timeevolution 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, SMstable 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. Twolevel 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 selfadjoint elliptic operator of second order is considered in Section 2. Using finite element approximations in space, we formulate a Cauchy problem for the corresponding differentialoperator equation. In Section 3, we consider standard twolevel 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 initialboundary value parabolic problem.
2 Problem statement
Let be a bounded domain () with a piecewise smooth boundary . Define an elliptic operator such that
(1) 
The operator is defined on the set of functions that satisfy on the boundary the following conditions:
(2) 
Coefficients , and are smooth functions in and
We consider the Cauchy problem
(3) 
(4) 
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
(5) 
Here is the following bilinear form:
where
with a constant . In view of (4), we put
(6) 
The variational (weak) formulation of the problem (1)–(4) consists in finding , that satisfies (5), (6) with the condition (2) on the boundary.
For the solution of the problem (5), (6
), we have the a priori estimate
(7) 
To show this, we put in (5) and get
Taking into account the lower bound for the bilinear form , we get
From this inequality, in view of (6), it follows that the estimate (7) holds.
For numerical solving the initialboundary 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
(8) 
where is the identity operator.
For the problem (3), (4), we put into the correspondence the operator equation for :
(9) 
(10) 
where with denoting projection onto . Similarly to (7), we establish the estimate
(11) 
To solve the problem (9), (10), the separation of variables method is applied. For the spectral problem
we have real eigenvalues
The eigenfunctions form a basis in . Therefore, for , we have
In view of this, the solution of the problem (9), (10) is represented in the form
(12) 
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 SMstable 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 twolevel 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 twolevel scheme:
(13) 
with the initial conditions
(14) 
where is a parameter (weight). If , then (13), (14) becomes the explicit scheme, for , we obtain the fully implicit scheme, whereas corresponds to the symmetric (Crank–Nicolson) scheme.
Similarly to (12), the solution of the problem (13), (14) can be represented as
(15) 
Substitution of (15) into (13) results in the following relation for an individual mode:
where
The scheme (13), (14) belongs to the class of FMES, if the following condition holds:
(16) 
For any , the equation
does not hold. Thus, for the schemes with weights (13), (14), there is not any fundamental mode exact scheme.
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 timevariation of the fundamental mode.
In our previous work afanas2013unconditionally , special unconditionally stable schemes were constructed for convectiondiffusion 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.
For the solution of the problem (9), (10), we apply the representation
(17) 
Substitution of (17) into (9), (10) results in the problem
(18) 
(19) 
where
(20) 
and
For the problem (18), (19), we use the standard scheme with weights:
(21) 
(22) 
In view of (17), suppose that . Thus, we arrive at the scheme
(23) 
with the initial condition
(24) 
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:
(25) 
Using the representation (15), we obtain
Thus, the condition (16) is fulfilled and this scheme is of FMEStype.
The nonstandard scheme (20), (23), (24) is unconditionally stable under the standard restriction on the weight . The proof is based on the wellknown (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
Using the relation between and , we obtain for the scheme (23), (24) the a priori estimate
(26) 
This estimate is consistent with the estimate (11) for the problem (9), (10), since . Thus, we arrive at the following statement.
4 Schemes based on Padé approximations
Twolevel finite difference schemes of higher approximation order for timedependent 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 RungeKutta method HairerWanner2010 ; Butcher2008 ; DekkerVerwer1984 .
The solution of the problem (9), (10) can be represented as
and so
For the corresponding difference scheme, we have
(27) 
with the operator of transition to a new time level . For (27), we have the representation (15), where now 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
(28) 
where and are polynomials of degrees and , respectively. These polynomials have BakerGravesMorris1996 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
The fundamental mode exact scheme based on Padé approximations can be represented in the form (27), where in view of (20) and (28), we have
(29) 
In fact, this means that the transformation (25) is realized. Thus, we can formulate the following statement.
Theorem 2
Let us highlight the most interesting variants of the constructed schemes of FMEStype. The scheme (27), (29) with corresponds to the fully implicit scheme ( in(23), (24)). The scheme (27), (29) with is associated with the CrankNicolson scheme ( in (23), (24)).
Only the schemes with , i.e.,
are SMstable Vabischevich2010b . In this case, the twolevel difference scheme for the problem (9), (10) is
(30) 
where the function
is a truncated Taylor series for .
In the class of schemes (30), in addition to the fully implicit scheme of the firstorder approximation in time (), we highlight the scheme with . In this case, the approximate solution is determined from
(31) 
It has, like the CrankNicolson scheme ( in (23), (24)), the secondorder approximation in time, but, as an SMstable scheme, it reproduces more accurately the behavior of higher modes.
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. Neumanntype or Robintype 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. Piecewiselinear 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 welldesigned algorithms and relevant opensource 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 KrylovSchur 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
(32) 
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 .
1  5.48146728860  5.41108707044  5.40129163974 

2  4.62435922303  4.54593036350  4.54096438823 
3  4.61225694503  4.53300663413  4.52815998124 
4  4.61203189452  4.53275210569  4.52790926702 
5  4.61202756655  4.53274692872  4.52790419871 
6  4.61202748265  4.53274682270  4.52790409555 
7  4.61202748102  4.53274682052  4.52790409345 
8  4.61202748099  4.53274682048  4.52790409341 
9  4.61202748099  4.53274682048  4.52790409340 
10  4.61202748099  4.53274682048  4.52790409340 

The timedependent 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 FMEStype 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.
Acknowledgements
The publication was financially supported by the Ministry of Education and Science of the Russian Federation (the Agreement # 02.a03.21.0008).
References
 (1) W. H. Hundsdorfer, J. G. Verwer, Numerical Solution of Timedependent Advectiondiffusionreaction 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: Steadystate and Timedependent Problems, Society for Industrial Mathematics, 2007.
 (5) E. Hairer, G. Wanner, Solving Ordinary Differential Equations II: Stiff and Differentialalgebraic 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 RungeKutta Methods for Stiff Nonlinear Differential Equations, NorthHolland, 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 reactiondiffusion equations, Numerical Methods for Partial Differential Equations 15 (2) (1999) 201–214.
 (15) R. Anguelov, P. Kama, J.S. Lubuma, On nonstandard 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. HernandezMartinez, H. Puebla, F. ValdesParada, J. AlvarezRamirez, 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, Twolevel finite difference scheme of improved accuracy order for timedependent 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 convectiondiffusion 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. GravesMorris, 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.