Fundamental mode exact schemes for unsteady problems

by   Petr N. Vabishchevich, et al.

The problem of increasing the accuracy of an approximate solution is considered for boundary value problems for parabolic equations. For ordinary differential equations (ODEs), nonstandard finite difference schemes are in common use for this problem. They are based on a modification of standard discretizations of time derivatives and, in some cases, allow to obtain the exact solution of problems. For multidimensional problems, we can consider the problem of increasing the accuracy only for the most important components of the approximate solution. In the present work, new unconditionally stable schemes for parabolic problems are constructed, which are exact for the fundamental mode. Such two-level schemes are designed via a modification of standard schemes with weights using Padé approximations. Numerical results obtained for a model problem demonstrate advantages of the proposed fundamental mode exact schemes.


Numerical solution of the Cauchy problem for Volterra integrodifferential equations with difference kernels

We consider the problems of the numerical solution of the Cauchy problem...

A Comparative Study of Iterative Riemann Solvers for the Shallow Water and Euler Equations

The Riemann problem for first-order hyperbolic systems of partial differ...

Approximate solution of the Cauchy problem for a first-order integrodifferential equation with solution derivative memory

We consider the Cauchy problem for a first-order evolution equation with...

Comparative study of Three Numerical Schemes for Fractional Integro differential Equations

This paper presents a comparative study three numerical schemes such as ...

Exact splitting methods for kinetic and Schrödinger equations

In [8], some exact splittings are proposed for inhomogeneous quadratic d...

Random forward models and log-likelihoods in Bayesian inverse problems

We consider the use of randomised forward models and log-likelihoods wit...

Accuracy of Difference Schemes in Electromagnetic Applications: a Trefftz Analysis

The paper examines local approximation errors of finite difference schem...

1 Introduction

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 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 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


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


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 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.

For the problem (3), (4), we put into the correspondence the operator equation for :


where with denoting -projection onto . Similarly to (7), we establish the estimate


for the solution of the problem (9), (10).

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


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


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


Substitution of (15) into (13) results in the following relation for an individual mode:


The scheme (13), (14) belongs to the class of FMES, if the following condition holds:


For any , the equation

does not hold. Thus, for the schemes with weights (13), (14), there is not any fundamental mode exact scheme.

Figure 1: Weight value for FMES.

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.

For the solution of the problem (9), (10), we apply the representation


Substitution of (17) into (9), (10) results in the problem





For the problem (18), (19), we use the standard scheme with weights:


In view of (17), suppose that . Thus, we arrive at the scheme


with the initial condition


This is a new scheme for the initial problem (9), (10).

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

Using the relation between and , we obtain for the scheme (23), (24) the a priori estimate


This estimate is consistent with the estimate (11) for the problem (9), (10), since . Thus, we arrive at the following statement.

Theorem 1

The difference scheme (20), (23), (24) is a fundamental mode exact scheme. It is unconditionally stable for and the estimate (26) holds.

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 .

The solution of the problem (9), (10) can be represented as

and so

For the corresponding difference scheme, we have


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


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

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


In fact, this means that the transformation (25) is realized. Thus, we can formulate the following statement.

Theorem 2

The difference schemes (20), (24), (27), (29) based on Padé approximations belong to the class of fundamental mode exact schemes.

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.,

are SM-stable Vabischevich2010b . In this case, the two-level difference scheme for the problem (9), (10) is


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


It has, like the Crank-Nicolson scheme ( in (23), (24)), the second-order approximation in time, but, as an SM-stable scheme, it reproduces more accurately the behavior of higher modes.

5 Numerical experiments

Figure 2: Computational domain, coefficients of the equation and the boundary conditions

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, 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 .

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

Table 1: Iterative approximation to the first eigenvalue
Figure 3: Fundamental eigenfunction.

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.

Figure 4: Time-variation the amplitude of the fundamental mode.

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.

Figure 5: Accuracy of the standard fully implicit scheme.
Figure 6: Accuracy of the nonstandard fully implicit scheme.
Figure 7: Accuracy of the standard fully implicit scheme for various values of .
Figure 8: Accuracy of the nonstandard fully implicit scheme for various values of .


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.