On asymptotic global error estimation and control of finite difference solutions for semilinear parabolic equations

Abstract

The aim of this paper is to extend the global error estimation and control addressed in Lang and Verwer [SIAM J. Sci. Comput. 29, 2007] for initial value problems to finite difference solutions of semilinear parabolic partial differential equations. The approach presented there is combined with an estimation of the PDE spatial truncation error by Richardson extrapolation to estimate the overall error in the computed solution. Approximations of the error transport equations for spatial and temporal global errors are derived by using asymptotic estimates that neglect higher order error terms for sufficiently small step sizes in space and time. Asymptotic control in a discrete $L_2$-norm is achieved through tolerance proportionality and uniform or adaptive mesh refinement. Numerical examples are used to illustrate the reliability of the estimation and control strategies.

Introduction

We consider semilinear parabolic partial differential equations

$${\partial }_{t}u(t,x)=L(t,x)u(t,x)+g(t,x,u(t,x))\text{,}t\in (0,T],x\in \Omega \subset {\mathbb{R}}^d\text{,}$$

in $d\in \mathbb{N}$ space dimensions, where $L$ is an elliptic operator, and assume that an appropriate system of boundary conditions and the initial condition

$$u(0,x)=u_0\left(x\right)\text{,}x\in \overline{\Omega }$$

are given. The initial boundary value problem is assumed to be well posed and to have a unique continuous solution $u(t,x)$.

The method of lines is used to solve (1) numerically. We first discretize the PDE in space by means of finite differences of order $q>1$ on a (possibly non-uniform) spatial mesh ${\Omega }_h$ and solve the resulting system of ODEs using existing time integrators. For simplicity, we shall assume that this system of time-dependent ODEs can be written in the general form

$$U_h^{\prime }\left(t\right)=F_h(t,U_h\left(t\right))\text{,}t\in (0,T]\text{,}$$$$U_h\left(0\right)=U_{h,0}\text{,}$$

with a unique solution vector $U_h\left(t\right)$ being a grid function on ${\Omega }_h$. Let

$$R_h:u(t,\cdot )\to \left(R_{h}u\right)\left(t\right)$$

be the usual restriction operator defined by $\left(R_{h}u\right)\left(t\right)={(u(t,x_1),\dots ,u(t,x_N))}^T$, where $x_i\in {\Omega }_h$ and $N$ is the number of all mesh points. Then we take as initial condition $U_{h,0}=R_{h}u\left(0\right)$.

To simplify the following derivations, we assume that $F_h$ is given by

$$F_h(t,U_h)=L_h\left(t\right)U_h+G_h(t,U_h)$$

with a finite difference approximation $L_h$ of $L$, and $G_h(t,R_{h}u)=R_{h}g(t,\cdot ,u(t,\cdot ))$.

To solve the initial value problem (3), we apply a numerical integration method of order $p\ge 1$ at a certain time grid

$$0=t_0<t_1<\cdots <t_n<\cdots <t_{M-1}<t_M=T\text{,}$$

using local control of accuracy. This yields approximations $V_h\left(t_n\right)$ to $U_h\left(t_n\right)$, which may be calculated for other values of $t$ by using a suitable interpolation method provided by the integrator. The global time error is then defined by

$$e_h\left(t\right)=V_h\left(t\right)-U_h\left(t\right)\text{.}$$

Numerical experiments in  [1] for ODE systems have shown that classical global error estimation based on the first variational equation is remarkably reliable. In addition, having the property of tolerance proportionality, that is, there exists a linear relationship between the global time error and the local accuracy tolerance, $e_h\left(t\right)$ can be successfully controlled by a second run with an adjusted local tolerance. Numerous techniques to estimate global errors are described in  [2]. A comparison of various adaptive grid methods for partial differential equations and implementation issues are presented in  [3], [4].

In order for the method of lines to be used efficiently, it is necessary to take also into account the spatial discretization error. Defining the spatial discretization error by

$${\eta }_h\left(t\right)=U_h\left(t\right)-\left(R_{h}u\right)\left(t\right)\text{,}$$

the vector of overall global errors $E_h\left(t\right)=V_h\left(t\right)-\left(R_{h}u\right)\left(t\right)$ may be written as sum of the global time and spatial error, that is,

$$E_h\left(t\right)=e_h\left(t\right)+{\eta }_h\left(t\right)\text{.}$$

We assume that $u(t,x)$ is $(q+2)$-times continuously differentiable with respect to $x$ and $(p+1)$-times continuously differentiable with respect to $t$. Then, with maximum step sizes $h_{\max }$ in space and ${\tau }_{\max }={max}_{i=0,\dots ,M-1}(t_{n+1}-t_n)$ in time it holds for the global space and time error that $\Vert {\eta }_h\left(t_n\right)\Vert =\mathcal{O}\left(h_{\max }^q\right)$ and $\Vert e_h\left(t_n\right)\Vert =\mathcal{O}\left({\tau }_{\max }^p\right),n=1,\dots ,M$, respectively.

Although a posteriori error estimates and adaptive algorithms for the efficient solution of parabolic problems are well established (see e.g.  [5], [6] and references therein), the separation of global time and spatial discretization errors is still a challenge. First experiences to estimate and balance the spatial discretization error and the error due to time integration of the ODEs within the method of lines have been made by Schönauer, Schnepf, and Raith  [7]. In their control strategy, the spatial mesh is initially chosen and remains fixed. The spatial truncation error is designated to be the level to which the local time error must be adapted. Lawson, Berzins, and Dew  [8] proposed to additionally control the local time error with respect to the contribution of the existing error from the previous time steps to the global error at the end of the next time step. The error in time is enabled to vary in relation to the spatial discretization error, ensuring that the method of lines with a fixed spatial mesh is being used efficiently. A successful attempt to assess and to equilibrate the individual discretization errors with respect to a given quantity of interest has been made by Schmich and Vexler  [9]. An adjoint linear parabolic problem has to be solved backwards in time to derive useful error bounds, which are used to enhance the resolution in time and space to meet a user-prescribed accuracy tolerance.

It is the purpose of this paper to present a new asymptotic error control strategy for the global errors $E_h\left(t\right)$, based on asymptotic estimates. We will mainly focus on reliability. So our aim is to provide error estimates ${\tilde{E}}_h\left(t\right)\approx E_h\left(t\right)$ which are not only asymptotically exact, but also work reliably for moderate tolerances, that is for relatively coarse discretizations. Approximations of the error transport equations for spatial and temporal global errors are derived by using asymptotic estimates that neglect higher order error terms for sufficiently small step sizes in space and time. The approximate global errors are measured in discrete $L_2$-norms. A priori bounds for the global error in such norms are well known, see e.g.  [10], [11]. However, reliable a posteriori error estimation and efficient control of the accuracy of the solution numerically computed to an imposed tolerance level are still challenging. We achieve asymptotic global error control by iteratively improving the temporal and spatial discretizations according to asymptotic estimates of $e_h\left(t\right)$ and ${\eta }_h\left(t\right)$. The global time error is estimated and controlled along the way fully described in  [1]. To estimate the global spatial error, we follow an approach proposed in  [12] (see also  [8]) and use Richardson extrapolation to set up a linearized error transport equation. Both strategies have to be combined in the right manner in order to make sure that they work reliably. Therefore, we have developed an appropriate control rule for the global spatial error. To control the overall global error more efficiently, we also consider a new fully space–time adaptive approach.

Throughout the paper we will use the terms ‘approximation’ and ‘estimation’ in the sense of asymptotic estimates, i.e., estimates that involve the Landau symbol $\mathcal{O}$.

The outline of this paper is as follows: In Section  2, we will linearize the transport equations for the global spatial and the global time error. These contain the residual time error and the spatial truncation error, which are approximated in Sections  3 Approximation of the residual time error, 4 Approximation of the spatial truncation error. In Section  5 we describe the discretization formulas used to approximate the solutions of the error transport equations, as well as the strategies used to adaptively adjust the time step size and the spatial mesh in dependence on the residual time error and the spatial truncation error. Now that we have approximations to the global time and global spatial error, Section  6 suggests strategies to adapt the local tolerances such that in further runs first the global time error and then the global spatial error respect some global tolerances provided by the user. Finally, numerical examples and a summary are given in Sections  7 Numerical illustrations, 8 Summary.