Arbitrary high order A-stable and B-convergent numerical methods for ODEs via deferred correction

TL;DR

This paper develops a sequence of high-order, A-stable, B-convergent numerical methods for solving first-order ODEs using recursive deferred correction schemes based on the implicit midpoint method.

Contribution

It introduces a new family of high-order deferred correction schemes that are A-stable and B-convergent, with proven order enhancement through correction steps.

Findings

Numerical experiments confirm high accuracy of schemes DC2 to DC10.

Theoretical order of accuracy is achieved in practice.

Schemes demonstrate satisfactory stability on stiff and non-stiff ODEs.

Abstract

This paper presents a sequence of deferred correction (DC) schemes built recursively from the implicit midpoint scheme for the numerical solution of general first order ordinary differential equations (ODEs). It is proven that each scheme is A-stable, satisfies a B-convergence property, and that the correction on a scheme DC(2j) of order 2j of accuracy leads to a scheme DC2j+2 of order 2j+2. The order of accuracy is guaranteed by a deferred correction condition. Numerical experiments with standard stiff and non-stiff ODEs are performed with the DC2, ..., DC10 schemes. The results show a high accuracy of the method. The theoretical orders of accuracy are achieved together with a satisfactory stability.

Tables

Table 1. Table 1: Coefficients of the approximations ( 13 )-( 14 ) for j = 1 , 2 , 3 , 4 𝑗 1 2 3 4 j=1,2,3,4

$j$	$c_2^j$	$c_3^j$	$c_4^j$	$c_5^j$	$c_6^j$	$c_7^j$	$c_8^j$	$c_9^j$
1	$\frac{9}{8}$	$\frac{9}{8}$
2	$\frac{25}{8}$	$\frac{125}{24}$	$\frac{125}{128}$	$\frac{125}{128}$
3	$\frac{49}{8}$	$\frac{343}{24}$	$\frac{637}{128}$	$\frac{13377}{1920}$	$\frac{1029}{1024}$	$\frac{1029}{1024}$
4	$\frac{81}{8}$	$\frac{243}{8}$	$\frac{1917}{128}$	$\frac{17253}{640}$	$\frac{7173}{1024}$	$\frac{64557}{7168}$	$\frac{32733}{32768}$	$\frac{32733}{32768}$

Table 2. Table 2: Absolute error (order of convergence) for the Bernoulli problem.

$k$	DC2	DC4	DC6	DC8	DC10
1	0.18	1.7e-2	1.8e-4	2.3e-4	1.3e-4
2.03e-3	3.71e-2 (0.26)	6.16e-4 (0.53)	7.14e-5 (0.14)	1.47e-6 (0.81)	9.42e-7 (0.79)
1.00e-4	1.92e-3 (0.98)	2.93e-5 (1.01)	4.31e-6 (0.93)	3.72e-7 (0.45)	5.78e-8 (0.94)
1.00e-5	2.22e-5 (1.94)	1.30e-7 (2.35)	3.92e-9 (3.04)	1.9e-10 (3.27)	1.1e-11 (3.73)
5.00e-6	5.55e-6 (2.0)	1.04e-8 (3.70)	1.4e-10 (3.70)	4.4e-12 (5.50)	4.4e-13 (4.64)
3.33e-6	2.46e-6 (1.99)	2.59e-9 (3.33)	1.6e-11 (5.31)	4.5e-13 (5.63)	2.0e-13 (2.02)
2.25e-6	1.39e-6 (1.99)	8.7e-10 (3.79)	3.3e-12 (5.54)	4.2e-13 (0.16)	4.2e-13 (-2.66)
$k$	BDF2	BDF4	BDF6	RK4
1	0.14	0.83	6.1e-2	–
2.03e-3	4.3e-2 (0.19)	2.5e-2 (0.19)	1.9e-3 (0.19)	–
1.00e-4	6.61e-3 (0.62)	2.98e-3 (0.71)	1.79e-3 (0.79)	1.27e-3
1.00e-5	2.59e-4 (1.41)	1.92e-5 (2.19)	3.15e-6 (2.76)	4.91e-8 (4.41)
5.00e-6	7.29e-5 (1.91)	1.92e-6 (3.58)	1.35e-7 (5.11)	2.53e-9 (4.28)

Table 3. Table 3: Absolute error (order of convergence) for the oscillatory problem.

$k$	DC2	DC4	DC6	DC8	DC10
5.00e-2	3418	456.26	42.665	3.2350	0.2132
2.50e-2	790.2 (2.1)	25.351 (4.2)	0.5959 (6.2)	1.17e-2 (8.1)	1.9e-4 (10.1)
1.25e-2	193.8 (2.0)	1.5493 (4.0)	9.17e-3 (6.0)	5.28e-5 (7.8)	2.79e-6 (6.1)
6.25e-3	48.23 (2.0)	9.67e-2 (4.0)	1.4e-4 (5.99)	2.78e-6 (0.0)	2.78e-6 (0.0)
1.56e-3	3.010 (2.0)	3.8e-4 (3.99)	4.72e-6 (2.5)	4.67e-6 (-0.3)	4.7e-6 (-0.3)
$k$	BDF2	BDF4	BDF6	rkf	stiff
1.56e-3	22026.46	14836.76	5578.40	22026.46	2636.00

Table 4. Table 4: Absolute error (order of convergence) for the first component of the solution for B ​ 5 𝐵 5 B5 modified

$k$	DC2	DC4	DC6	DC8	DC10
2.000e-5	0.2152	6.51e-2	2.22e-2	8.00e-3	2.98e-3
5.000e-6	1.35e-2 (2)	2.59e-4 (4)	5.59e-6 (6)	1.27e-7 (8)	2.97e-9 (10)
2.500e-6	3.38e-3 (2)	1.62e-5 (4)	8.74e-8 (6)	4.9e-10 (8)	2.9e-12 (10)
1.250e-6	8.47e-4 (2)	1.01e-6 (4)	1.36e-9 (6)	1.9e-12 (8)	7.4e-14 (5.3)
3.125e-7	5.29e-5 (2)	4.00e-9 (4)	3.6e-13 (6)	7e-14 (2.4)	6.3e-14
6.250e-8	2.11e-6 (2)	6.3e-12 (4)	6.02e-13	2.33e-13	1.19e-13
$k$	BDF2	BDF4	BDF6	rkf	stiff
1.25e-6	3.38e-3	7.94e-8	2.3e-12	2.36e-6	6.6e-10

Table 5. Table 5: Absolute error (order of convergence) for the problem E5

$k$	DC2	DC4	DC6	DC8	DC10
100	2.79e-07	5.34e-08	8.31e-09	4.26e-09	1.04e-09
	8.30e-12	9.68e-13	6.86e-14	6.14e-14	1.66e-14
	4.47e-13	5.31e-14	3.28e-15	3.40e-15	8.42e-16
	7.85e-12	9.14e-13	6.54e-14	5.81e-14	1.57e-14
50	7.52e-08(1.89)	1.02e-08(2.38)	1.56e-09(2.41)	8.53e-11(5.64)	4.92e-11(4.41)
	1.96e-12(2.08)	6.46e-14(3.90)	3.16e-14(1.12)	2.94e-15(4.38)	5.07e-16(5.03)
	1.07e-13(2.06)	3.73e-15(3.83)	1.61e-15(1.02)	2.21e-16(3.94)	9.78e-17(3.11)
	1.86e-12(2.08)	6.14e-14(3.89)	3.00e-14(1.12)	2.85D-15(4.35)	4.09D-16(5.26)
10	3.16e-09(1.99)	2.37e-11(4.03)	5.26e-13(5.23)	1.28e-14(6.72)	4.51e-16(8.89)
	7.77e-14(1.99)	2.79e-16(3.68)	3.02e-18(5.74)	1.15e-19(7.94)	7.28e-21(8.09)
	4.31e-15(1.97)	7.08e-17(1.79)	5.91e-17(0.24)	6.27e-17(0.12)	6.84e-17(0.09)
	7.34e-14(1.99)	3.20e-16(3.37)	6.18e-17(1.79)	6.28e-17(0.57)	6.84e-17(0.11)
$k$	BDF2	BDF4	BDF6	RK4	stiff
10	5.7e-8	6.6e-10	3.5e-11	2.03e-16	1.29e-16

Table 6. Table 6: Absolute error (order of convergence) for Robertson problem

$k$	DC2	DC4	DC6	DC8	DC10
0.5	3.63e-5	4.46e-6	2.08e-6	2.91e-6	3.09e-6
	3.63e-5	4.46e-6	2.08e-6	2.91e-6	3.09e-6
	7.12e-5	4.37e-7	1.02e-7	4.12e-7	4.26e-7
1/300	4.7e-9 (1.8)	1.09e-9 (1.7)	4.0e-10 (1.7)	3.0e-10 (1.9)	2.0e-10 (1.9)
	7.4e-9 (1.7)	2.23e-8 (1.1)	4.16e-8 (0.8)	2.9e-8 (0.9)	2.5e-8 (0.9)
	4.7e-9 (1.9)	2.12e-8 (0.6)	4.12e-8 (0.6)	2.8e-8 (0.5)	2.5e-8 (0.6)
1/600	1.0e-9 (2.2)	1.5e-10 (2.8)	1.0e-12 (8.6)	9.9e-13 (8.)	7.5e-13 (8.2)
	5e-13 (14.)	3.0e-14 (19.6)	2.0e-16 (27.7)	2.0e-16 (27.1)	3.0e-16 (26.1)
	1.0e-9 (2.2)	1.5e-10 (7.1)	1.0e-12 (15.3)	9.9e-13 (15)	4.0e-13 (15.8)
1/6000	9.24e-12	7.31e-14	1.48e-14	4.57e-14	–
	5.38e-15	0.	0.	0.	–
	9.25e-12	2.07e-13	1.36e-13	8.27e-14	–
$k$	BDF2	BDF4	BDF6	RK4	stiff
0.5	5.3e-4	3.6e-5	4.1e-6	–	7.76e-13
1/600	2.8e-6	1.2e-6	6.9e-7	–	7.28e-13

Table 7. Table 7: Absolute error (order of convergence) for the van der Pol’s equation

$k$	DC2	DC4	DC6	DC8	DC10
3.75e-5	3.0089	2.9999	2.9440	0.1838	3.12e-3
	1322.9	1327.5	1320.6	197.79	3.26792
1.50e-5	2.9769 (0)	2.9999 (0)	0.1080 (3.6)	1.90e-4 (7.5)	5.1e-5 (4.5)
	1333.3 (0)	1330.3 (0)	113.69 (2.7)	0.18281 (7.6)	5.1e-2 (4.5)
7.50e-6	2.8706 (0)	2.6947 (0)	1.60e-3 (6.0)	1.74e-6 (6.7)	1.27e-5 (1.9)
	1327.4 (0)	1286.5 (0)	1.6349 (6.1)	1.80e-3 (6.7)	1.29e-2 (1.9)
1.875e-6	0.74(0.9)	0.339 (1.5)	2.50e-7 (6.3)	–	2.88e-7 (2.7)
	659. (0.5)	373.2 (0.9)	2.91e-4 (6.2)	–	2.92e-4 (2.7)
–	stiff	rkf
	2.16e-6	3.54e-2
	3.48e-3	64.76

Equations

$$\left\{\begin{array}[]{cccc}\displaystyle\frac{du}{dt}&=&F(t,u),&~{}~{}t\in[0,T],\\ u(0)&=&u_{0},&\end{array}\right.$$

$$Du(t_{n+1/2})=\frac{u(t_{n+1})-u(t_n)}{k},$$

$$D_{+}u(t_n)=\frac{u(t_{n+1})-u(t_n)}{k},$$

$$D_{-}u(t_n)=\frac{u(t_n)-u(t_{n-1})}{k},n\ge 1.$$

$$Eu(t_{n+1/2})=\hat{u}(t_{n+1})=\frac{u(t_{n+1})+u(t_n)}{2}.$$

$$(D_{+}{D}_{-}{)}^{m}u(t_n)=k^{-2m}\sum _{i=0}^{2m}(-1{)}^i\left(\genfrac{}{}{0ex}{}{2m}{i}\right)u(t_{n+m-i}),$$

$$D_{-}(D_{+}{D}_{-}{)}^{m}u(t_n)=k^{-2m-1}\sum _{i=0}^{2m+1}(-1{)}^i\left(\genfrac{}{}{0ex}{}{2m+1}{i}\right)u(t_{n+m-i}),$$

$$\Vert D_{+}^{m_1}{D}_{-}^{m_2}u(t_n)\Vert \le \underset{0\le t\le T}{\max }\Vert \frac{d^{m_1+m_2}u}{d{t}^{m_1+m_2}}(t)\Vert ,$$

$$D{u}^{n+1/2}=D_{+}{u}^n=D_{-}{u}^{n+1}=\frac{u^{n+1}-u^n}{k},$$

$$E{u}^{n+1/2}={\hat{u}}^{n+1}=\frac{u^{n+1}+u^n}{2}.$$

$$\frac{du}{dt}(t_{n+1/2})=\frac{u(t_{n+1})-u(t_n)}{k}-\sum _{i=1}^j{c}_{2i+1}{k}^{2i}(D_{+}{D}_{-}{)}^{i}Du(t_{n+1/2}))+O(k^{2j+2})$$

$$u(t_{n+1/2})=\frac{u(t_{n+1})+u(t_n)}{2}-\sum _{i=1}^j{c}_{2i}{k}^{2i}(D_{+}{D}_{-}{)}^{i}Eu(t_{n+1/2})+O(k^{2j+2}),$$

$$\frac{u^{n+1}-u^n}{k}-$$

$$=F\left(t_{n+1/2},\frac{u^{n+1}+u^n}{2}-\sum _{i=1}^j{c}_{2i}{k}^{2i}(D_{+}{D}_{-}{)}^{i}E{u}^{n+1/2}\right).$$

$$\frac{u^{2,n+1}-u^{2,n}}{k}=F\left(t_{n+1/2},\frac{u^{2,n+1}+u^{2,n}}{2}\right),u^{2,0}=u_0.$$

$$\frac{u^{2j+2,n+1}-u^{2j+2,n}}{k}-\sum _{i=1}^j{c}_{2i+1}{k}^{2i}(D_{+}{D}_{-}{)}^{i}D{u}^{2j,n+1/2}$$

$$=F\left(t_{n+1/2},\frac{u^{2j+2,n+1}+u^{2j+2,n}}{2}-\sum _{i=1}^j{c}_{2i}{k}^{2i}(D_{+}{D}_{-}{)}^{i}E{u}^{2j,n+1/2}\right),$$

$$u^{2j+2,0}=u_0.$$

$$\frac{u^{2j+2,n+1}-u^{2j+2,n}}{k}-k^{-1}\sum _{i=1}^j{c}_{2i+1}^j{k}_j^{2i+1}(D_{+}{D}_{-}{)}^{i}D{\bar{u}}^{2j,(2j+1)n+j+1/2}$$

$$=F\left(t_{n+1/2},E{u}^{2j+2,n+1/2}-\sum _{i=1}^j{c}_{2i}^j{k}_j^{2i}(D_{+}{D}_{-}{)}^{i}E{\bar{u}}^{2j,(2j+1)n+j+1/2}\right),$$

$$u^{2j+2,0}=u_0.$$

$$u^{\prime }(t_{n+1/2})=\frac{u(t_{n+1})-u(t_n)}{k}-\frac{1}{k}\sum _{i=1}^j{c}_{2i+1}^j{k}_j^{2i+1}D(D_{+}{D}_{-}{)}^{i}u({\tau }_{j+1/2})+O(k_j^{2j+2})$$

$$u(t_{n+1/2})=\frac{u(t_{n+1})+u(t_n)}{2}-\sum _{i=1}^j{c}_{2i}^j{k}_j^{2i}(D_{+}{D}_{-}{)}^{i}Eu({\tau }_{j+1/2})+O(k_j^{2j+2}),$$

$$x-a_n^j-kF(t_{n+1/2},0.5x+b_n^j)=0,$$

$$\Vert (D_{+}{D}_{-})D(u^{2j,n+1/2}-u(t_{n+1/2}))\Vert +\Vert D_{+}{D}_{-}(u^{2j,n+1}-u(t_{n+1}))\Vert \le C{k}^{2j},$$

$$\Vert \sum _{i=1}^j{c}_{2i}{k}^{2i}{\left(D_{+}{D}_{-}\right)}^i\left(u^{2j,n}-u(t_n)\right)\Vert \le C{k}^{2j+2},$$

$$\Vert \sum _{i=1}^j(c_{2i+1}-c_{2i})k^{2i}{\left(D_{+}{D}_{-}\right)}^{i}D\left(u^{2j,n+1/2}-u(t_{n+1/2})\right)\Vert \le C{k}^{2j+2},$$

$$k^{2i}{\left(D_{+}{D}_{-}\right)}^i\left(u^{2j,n}-u(t_n)\right)=k^2\sum _{l=0}^{i-1}(-1{)}^l(\genfrac{}{}{0ex}{}{2i-2}{l})D_{+}{D}_{-}\left(u^{2j,n}-u(t_n)\right)$$

$$\Vert u^{2j+2,n}-u(t_n)\Vert \le C{k}^{2j+2},\text{mbox}forn=1,2,...,j,$$

$$\Vert F(t,x)-F(t,y)\Vert \le \mu \Vert x-y\Vert ,\forall (t,x,y)\in [0,T]\times X\times X.$$

$$\Vert u^{2j+2,n}-u(t_n)\Vert \le M,\text{mbox}foreachn=0,1,...,N.$$

$$\left(F(t,x)-F(t,y),x-y\right)\le \beta \Vert x-y{\Vert }^2,\forall (t,x,y)\in [0,T]\times X\times X.$$

Full text

††journal: BIT

1 Introduction

In MR2058857 ; kress2002deferred , Gustafsson and Kress introduced a new version of deferred correction (DC) strategy for the numerical solution of linear systems of ordinary differential equations (ODE) MR2058857 and initial boundary value problems kress2002deferred , under a monotonicity condition. Numerical experiments with one-dimensional linear parabolic and hyperbolic equations were performed and showed that the method is effective (orders 2, 4 and 6 of accuracy are achieved). We propose to extend the method from MR2058857 ; kress2002deferred to the time-discretization of more general time-evolution partial differential equations (PDEs). In this paper, we restrict to the case of the initial value problem (IVP)

[TABLE]

where the unknown $u$ is from $[0,T]$ into a Banach space $X$, $u_{0}$ is a given data and $F$ is a sufficiently differentiable function such that $u$ exists and is sufficiently differentiable. The main objective is to show the properties of the numerical method (consistency, stability, convergence and order of accuracy). A complete analysis of the DC method applied to reaction-diffusion equations leads to an arbitrary high order and unconditionally stable method (see koyaguerebo2020unconditionally ).

The DC method is used to improve the order of accuracy of numerical methods of lower order. This method is explored by many authors, e.g. schild1990gaussian ; auzinger2016encyclopedia ; MR2058857 ; kushnir2012highly ; hansen2011order ; dutt2000spectral ; daniel1967interated ; IntegralDC2010 . The method in daniel1967interated is an application of iterative deferred correction (IDC). The authors proved that an asymptotic improvement of order $p$ can be accomplished, from a scheme of order $p$, at each step of the IDC procedure, provided suitable finite difference operators are employed. Numerical experiments are performed with the IDC applied to the trapezoidal rule, Taylor-2 and Adams-Bashforth of order 2. The results are promising even though they point out some difficulties of the proposed algorithms: inaccuracy for “large” time step and no asymptotic improvement for high levels of correction. The approaches in kushnir2012highly ; hansen2011order ; dutt2000spectral ; auzinger2016encyclopedia ; MR2058857 ; IntegralDC2010 are quite similar and consist in a linear perturbation of a low order scheme. However, solving stiff problems (problems extremely hard to solve by standard explicit methods spijker1996stiffness ) is a challenge unfavorable for these methods. In particular, the method in kushnir2012highly , concerning a highly accurate solver for stiff ODEs, requires sufficiently small time steps for moderately stiff problems while convergence is reduced to order 2 for “very stiff” problems.

Our schemes are based on nonlinear perturbations (corrections) of the implicit midpoint rule and inherit the A-stable property of the trapezoidal rule MR0170477 at any stage of the correction. Starting from an approximation $\left\{u^{2,n}\right\}_{n=0}^{N}$ of the exact solution $u$ by the implicit midpoint rule on a uniform partition $0=t_{0}<t_{1}<\cdots<t_{N}=T$ of $[0,T]$, at the stage $j=1,2,\cdots$ of the correction we obtain an approximation $\left\{u^{2j+2,n}\right\}_{n=0}^{N}$ of $u$, expected to be of order $2j+2$ of accuracy, on the same partition. Each approximate solution $\left\{u^{2j,n}\right\}_{n=0}^{N}$ to be corrected is subject to a deferred correction condition (DCC) which guarantees the improvement of the order of accuracy. We prove that if $\left\{u^{2j,n}\right\}_{n=0}^{N}$ satisfies the DCC and its correction $\left\{u^{2j+2,n}\right\}_{n=0}^{N}$ converges to $u$ at the discrete points $0=t_{0}<t_{1}<\cdots<t_{N}=T$ (or is simply bounded, when $X$ is finite dimensional) then $\left\{u^{2j+2,n}\right\}_{n=0}^{N}$ approximates $u$ with order $2j+2$. Moreover, provided the function $F$ is Lipschitz with respect to its second variable or satisfies a one-sided Lipschitz condition, each $\left\{u^{2j,n}\right\}_{n=0}^{N}$ satisfies the DCC and then converges with order $2j$ of accuracy, for arbitrary positive integer $j$. We also prove that each DC scheme involving $\left\{u^{2j,n}\right\}_{n=0}^{N}$ is $B$-stable. The theory is illustrated by numerical tests, for the schemes of order 2, 4, …, 10.

The paper is organized as follows: in section 2 we recall some basic results from finite difference approximations and present the DC schemes; section 3 deals with the consistency of the method; the analysis of convergence and order of accuracy together with a B-convergence result are given in section 4; absolute stability is proved is section 5, and section 6 is devoted to numerical experiments.

2 Deferred correction schemes for the implicit midpoint rule

We suppose that $\displaystyle F\in C^{2p+2}\left([0,T]\times X,X\right)$, for a positive integer $p$, so that (1) has a unique solution $\displaystyle u\in C^{2p+3}\left([0,T],X\right)$. We simply denote by $\|\cdot\|$, the norm in the Banach space $X$. For a time step $k>0$, we denote $t_{n}=nk$ and $t_{n+1/2}=(n+1/2)k$, for each integer $n$. This implies that $t_{0}=0$. We consider the time steps $k$ such that $0=t_{0}<t_{1}<\cdots<t_{N}=T$ is a partition of $[0,T]$, for a non-negative integer $N$. The centered, forward and backward difference operators $D$, $D_{+}$ and $D_{-}$, respectively, related to $k$ and applied to $u$, are defined as follows:

[TABLE]

[TABLE]

and

[TABLE]

The average operator is denoted by $E$:

[TABLE]

The composition of $D_{+}$ and $D_{-}$ is defined recursively. They commute, that is $(D_{+}D_{-})u(t_{n})=(D_{-}D_{+})u(t_{n})=D_{-}D_{+}u(t_{n})$, and satisfy the identities

[TABLE]

and

[TABLE]

for each integer $m\geq 1$ such that $0\leq t_{n-m-1}\leq t_{n+m}\leq T$. We have the estimate

[TABLE]

provided $[t_{n-m_{2}},t_{n+m_{1}}]\subset[0,T]$ and $m_{1}+m_{2}\leq 2p+3$ (see (isaacson1966analysis, , p.249) or koyaguerebo2020finite ).

If $\left\{u^{n}\right\}_{n}$ is a sequence of approximation of $u$ at the discrete points $t_{n}$, the finite difference operators apply to $\left\{u^{n}\right\}_{n}$, and we define

[TABLE]

and

[TABLE]

From the centered finite difference approximation (see (koyaguerebo2020finite, , Thm 5) or hildebrand1974introduction ; chung2010computational ; dahlquist2008numerical ) we have

[TABLE]

and

[TABLE]

for each integer $j=1,2,\cdots,p$. These approximations lead to the schemes

[TABLE]

The schemes (7) are multi-steps and prone to stability restrictions. We resort to DC method to transform them into a sequence of one step schemes as follows: For $j=0$, we have the implicit midpoint rule

[TABLE]

For $j\geq 1$,

[TABLE]

[TABLE]

The scheme (LABEL:a27)-(10) has unknowns $u^{2j+2,n}$, $n=1,2,...,N$, and is deduced from (7) by substituting the unknown $u^{n}$ under the summation symbols by $u^{2j,n}$. The index $2j$ indicates that $\left\{u^{2j,n}\right\}_{n}$ is expected to be an approximation of the exact solution $u$ with order $2j$ of accuracy. We call the schemes (LABEL:a27)-(10) Deferred Correction of order $2j+2$ for the implicit midpoint rule, denoted DC2j+2.

Remark 1

The scheme (LABEL:a27)-(10), for $n=1,2,3,\cdots,j$, should involve unknowns $u^{2j,-1},...,u^{2j,-j}$ which represent approximate solutions of (1) at the discrete points $t=-k,...,-jk$, respectively. To avoid those approximations for $t<0$, we propose the following scheme which is efficient for the computation of $u^{2j+2,1},...,u^{2j+2,j}$, using only points within the solution interval $[0,T]$.

[TABLE]

[TABLE]

The finite difference operator in (11) are related to the time step $k_{j}=k/(2j+1)$. The approximations $\left\{\overline{u}^{2j,m}\right\}_{m}$ and $\left\{u^{2j,n}\right\}_{n}$ are computed from the same scheme, (8) or (LABEL:a27)-(10), but for the time steps $k_{j}$ and $k$, respectively. The scheme (11) results from the finite difference approximations

[TABLE]

and

[TABLE]

where $t_{n}=\tau_{0}<\tau_{1}<...<\tau_{2j+1}=t_{n+1}$, with $\tau_{m}=t_{n}+mk_{j}$, for $m=1,2,\cdots,2j+1$. Table 1 gives the coefficients ${c}^{j}_{i}$ for $j=1,2,3,4$.

Remark 2

Each $u^{2j+2,n+1}$, $n\geq j$, is an iterative solution of the system

[TABLE]

where $x$ is the unknown, and $a_{n}^{j}$ and $b_{n}^{j}$ are constants depending on $u^{2j+2,n}$ and $u^{2j,n+1+j},u^{2j,n+j},\cdots,u^{2j,n-j}$. The total number of vectors (in the solution space $X$) stored for the computation of $u^{2j+2,n+1}$ is $j^{2}+3j+1$: $u^{2j+2,n}$ and the $u^{2i,q}$, for $i=1,2,\cdots,j$, and $n+(j-i+1)(j+i)/2-2i\leq q\leq n+1+(j-i+1)(j+i)/2$.

Remark 3

From Remark 2, only the implicit midpoint rule, DC2, is an implicit Runge-Kutta (RK) methods. Starting with DC4, all the DC2j methods of the form (LABEL:a27)-(10) are not RK methods. For instance, $u^{4,n+1}$ depends on $u^{4,n}$ and some of the $u^{2,i}$, which $u^{2,i}$ evolve independently and are not stages computed from $u^{4,n}$. As we will see in Section 5, the analysis of A-stability, in particular the proof of lemma 3, shows that it is impossible to write a recurrence $u^{2j+2,n+1}=R(z)\,u^{2j+2,n}$ from (LABEL:a27) when $j\geq 1$, as one would get by applying any RK method to Dahlquist equation. This is the main ingredient behind the A-stability of our DC2j methods independently of the order of accuracy.

3 Deferred correction condition (DCC)

In this section we give a sufficient condition for the scheme (LABEL:a27)-(10) to achieve order $2j+2$ of accuracy. Hereafter, the letter $C$ will denote any constant independent from $k$, and that can be calculated explicitly in terms of known quantities. The exact value of $C$ may change. We have the following definition:

Definition 1

(Deferred Correction Condition) Let $u$ be the exact solution of the Cauchy problem (1). Given a positive integer $j$, a sequence $\left\{u^{2j,n}\right\}_{n=0}^{N}$ of approximations of $u$, at the discrete points $0=t_{0}<\cdots<t_{N}=T$, is said to satisfy the Deferred Correction Condition $(DCC)$ for the implicit midpoint rule if $\left\{u^{2j,n}\right\}_{n=0}^{N}$ approximates $u$ with order $2j$ of accuracy, and we have

[TABLE]

for $n=1,2,...,N-2$ and $k\leq k_{0}$, where $k_{0}>0$ is fixed and $C$ is a constant independent from $k$.

Remark 4

Condition (16) is equivalent to

[TABLE]

and

[TABLE]

for $n=j,j+1,\cdots,N-j$. This is due to the transform

[TABLE]

and a similar transform for $k^{i}\left(D_{+}D_{-}\right)^{i}D\left(u^{2j,n+1/2}-u(t_{n+1/2})\right)$.

We have the following result:

Theorem 3.1

Let $u$ be the exact solution of (1) and $\left\{u^{2j,n}\right\}_{n=0}^{N}$, $j=1,\dots,p$, a sequence of approximations of $u$ satisfying DCC for the implicit midpoint rule. Let $\left\{u^{2j+2,n}\right\}_{n=0}^{N}$ be the solution of (LABEL:a27)-(10) built from $\left\{u^{2j,n}\right\}_{n=0}^{N}$. We suppose that $u^{2j+2,1},...,u^{2j+2,j}$ are given and satisfy

[TABLE]

where $C$ is a constant independent from $k$. Furthermore, we suppose that one of the following four conditions holds:

$F$* is Lipschitz with respect to the second variable $x$: there exists $\mu\geq 0$ such that*

[TABLE]

$X$* is finite dimensional, and $\left\{u^{2j+2,n}\right\}_{n=0}^{N}$ remains close to $u$ in the sense that there exists $M>0$ such that*

[TABLE]

$X$* is infinite dimensional, and $\left\{u^{2j+2,n}\right\}_{n}$ converges to the exact solution $u$.*

$X$* is a Hilbert space with inner product $\left(.,.\right)$, and $F$ satisfies the following so-called one-sided Lipschitz condition, with a one-sided Lipschitz constant $\beta\in\mathbb{R}$:*

[TABLE]

Then $\left\{u^{2j+2,n}\right\}_{n}$ approximates $u$ with order $2j+2$ of accuracy, that is

[TABLE]

where $C$ is a constant depending only on $j$, $T$, DCC, a Lipschitz constant on $F$ and the derivatives of $u$ up to order $2j+3$, for time steps $k$ sufficiently small.

Proof

First we consider the case where the function $F=F(t,x)$ is Lipschitz with respect to the second variable $x$. Combining (1) and (LABEL:a27), we obtain the identity

[TABLE]

where $\Lambda^{j}$ and $\Gamma^{j}$ are finite difference operators defined for arbitrary integer $j\geq 1$ by

[TABLE]

and

[TABLE]

provided $u(t_{n\pm i})$ exists for $i=0,1,2,\cdots,j$. We have defined

[TABLE]

and

[TABLE]

From (5) we have

[TABLE]

and, since $F$ is differentiable and $u$ is sufficiently regular, we deduce from the mean value theorem and the approximation (6) that

[TABLE]

for each $n=0,1,\cdots,N$, where $C$ is a constant depending only on $j$, $T$, a Lipschitz constant from $F$ and the derivatives of $u$ up to order $2j+3$. The last two inequalities imply that

[TABLE]

Since the sequence $\left\{u^{2j,n}\right\}_{n}$ satisfies DCC, from Remark 4 we have

[TABLE]

From the Lipschitz condition on $F$ we have

[TABLE]

Substituting inequalities (26)-(28) in the identity (24), we deduce that

[TABLE]

and it follows from the triangle inequality that

[TABLE]

for $0\leq\mu k<2$. We then deduce by induction on $n$ that

[TABLE]

From hypothesis (19) and the DCC we have

[TABLE]

where $C$ is a constant independent from $k$. Moreover, the sequence $\left\{\left(\frac{2+\mu k}{2-\mu k}\right)^{n}\right\}_{n}$ is bounded above by $\exp(2\mu T/(2-\varepsilon))$, for $0\leq\mu k\leq\varepsilon<2$. Whence

[TABLE]

Finally, by the triangle inequality, identity (25) and DCC, we have

[TABLE]

where $C$ is a constant depending only on $j$, $T$, the DCC constant, $\mu$ and the derivatives of $u$ up to order $2j+3$.

Suppose that $\left\{u^{2j+2,n}\right\}_{n=0}^{N}$ satisfies (21) and $X$ is finite dimensional. We can write

[TABLE]

From (21) and the DCC there exists $k_{1}>0$ such that $0<k\leq k_{1}\leq k_{0}$ implies

[TABLE]

On the other hand, we have

[TABLE]

where

[TABLE]

It follows (28) for

[TABLE]

Since $F$ is differentiable and the set $\left\{x\in X:\|x\|\leq M+R_{j+1}+1\right\}$ is compact in the finite dimensional linear space $X$, the supremum exists and is finite. The theorem is then deduced from the case (i).

If $\left\{u^{2j+2,n}\right\}_{n}$ converges to the exact solution $u$, taking the DDC and the finite difference formula (6) into account, we have

[TABLE]

It follows from the continuity of $u\mapsto d_{u}F(t,u)$ that there exists $0<k_{2}\leq k_{0}$ such that $0<k\leq k_{2}$ implies

[TABLE]

The theorem, in this case, follows by taking $\mu=1+\max_{0\leq t\leq T}\|d_{u}F\left(t,u(t)\right)\|$ in (i).

Here we consider the case where $X$ is a Hilbert space and $F$ satisfies the monotonicity condition (22). Then, taking the inner product of the identity (24) with $\widehat{\Theta}^{2j+2,n+1}$, we deduce the inequality

[TABLE]

since, according to (22), we have

[TABLE]

Inequalities (26)-(27) together with the Cauchy-Schwartz inequality yield

[TABLE]

and

[TABLE]

where $C$ is a constant depending only on $j$, $T$, a Lipschitz constant on $F$ and the derivatives of $u$ up to order $2j+3$. Substituting the last three inequalities into (33), we obtain

[TABLE]

and we deduce from the identity

[TABLE]

and the inequality

[TABLE]

that

[TABLE]

The conclusion follows from the case (i), for $-2\leq\beta k<2$.

Remark 5

Theorem 3.1 shows that the correction may be applied for any other scheme satisfying DCC.

4 Convergence and order of accuracy

The aim of this section is to prove the following theorem:

Theorem 4.1

Let $u\in C^{2p+3}\left([0,T],X\right)$ be the exact solution of the problem (1). Suppose that one of the four conditions (i)-(iv) of Theorem 3.1 holds, with condition (ii) or (iii) holding for all $j=0,1,\cdots,p+1$. Then each sequence $\left\{u^{2j,n}\right\}_{n=0}^{N}$, $j=1,2,\cdots,p+1$, solution of the scheme (8) or (LABEL:a27)-(10), approximates $u$ with order $2j$ of accuracy. Furthermore, we have the estimate

[TABLE]

for $m=0,1,...,p-j$ and $n=m+j-1,m+j,...,N-j-m$, where $C$ is a constant depending only on $p$, $T$, and the derivatives of $u$ and $F$ up to order $2m+2j+1$ and $2m+2j-1$, respectively.

To prove this theorem we need Theorem 3.1 and the the following lemma:

Lemma 1

Let $\left\{u^{2,n}\right\}_{n=0}^{N}$ be the solution of the scheme (8). Suppose that one of the conditions (i), (iii) or (iv) of Theorem 3.1 holds, or $\left\{u^{2,n}\right\}_{n=0}^{N}$ is bounded in the sense of the condition (ii) of this theorem. Then $\left\{u^{2,n}\right\}_{n=0}^{N}$ approximates $u$ with order 2 of accuracy, and we have the inequality

[TABLE]

for $m=0,1,...,p$ and $n=m,m+1,...,N-m-1$, where $C$ is a constant depending only on $p$, $T$, and the derivatives of $u$ and $F$ up to order $2m+3$ and $2m+1$, respectively.

Proof (Proof of Lemma 1)

For the sake of simplification we suppose that $F=F(x)$. The general case can be handled by transforming (1) to an autonomous system. From the hypotheses of the Lemma, Theorem 3.1 implies that $\left\{u^{2,n}\right\}_{n=0}^{N}$ approximates $u$ with order two of accuracy:

[TABLE]

where $C$ is a constant depending only on $T$, $F$ and the derivatives of $u$ up to order 3. To establish (35) we proceed by induction on the integer $m=0,1,\cdots,p$.

Inequality (35) for $m=0$.

As in Theorem 3.1, we combine (1) and (8) and deduce the identity

[TABLE]

where

[TABLE]

and

[TABLE]

From Taylor’s formula with integral remainder and the estimate (4), there exists a function $g$ such that

[TABLE]

with

[TABLE]

for each nonnegative integers $m_{1}$ and $m_{2}$ such that $m_{1}+m_{2}\leq 2p$, where $C$ is a constant depending only on $T$, $F$, and the derivatives of $u$ up to order $m_{1}+m_{2}+3$. We can write

[TABLE]

where

[TABLE]

The last identities substituted into (37) yield

[TABLE]

Proceeding as in Theorem 3.1, we deduce from (36) and the regularity of $u$ that

[TABLE]

Therefore, taking the norm on both sides of (39), we deduce by the triangle inequality and the inequalities (36) and (38), for $m_{1}=m_{2}=0$, that

[TABLE]

where $C$ is a constant depending only on $T$ and the derivatives of $u$ and $F$ up to order 3 and 1, respectively. The last inequality combined with (36) implies that (35) holds for $m=0$.

Here we are going to prove that inequality (35) remains true for $m+1$, assuming that it holds for an arbitrary integer $m$ such that $0\leq m\leq p-1$.

We apply $\left(D_{+}D_{-}\right)^{m}D_{+}$ to (39) and obtain

[TABLE]

where we set

[TABLE]

The main difficulty is to bound $\left(D_{+}D_{-}\right)^{m}D_{+}h(t_{n+1})=D_{+}^{2m+1}h(t_{n+1-m})$. We have

[TABLE]

[TABLE]

where $d\tau^{i}=d\tau_{1}\cdots d\tau_{i}$, and

[TABLE]

It follows the general formula

[TABLE]

where $\alpha_{i}=(\alpha_{i}^{1},\cdots,\alpha_{i}^{i-1},\alpha_{i}^{i})\in\left\{1,2,\cdots,q\right\}^{i-1}\times\left\{0,1,\cdots,q-i+1\right\}$, and $L^{n,q}_{i,\alpha_{i}}$ is a linear combination, with properly chosen coefficients, of the quantities

[TABLE]

where $\beta_{i}=(\beta_{i}^{1},\cdots,\beta_{i}^{i-1},\beta_{i}^{i})\in\left\{1,2,\cdots,q\right\}^{i-1}\times\left\{0,1,\cdots,q-i+1\right\}$ with $\beta_{i}^{l}+\alpha_{i}^{l}\leq q-l+1$, for $l=1,\cdots,i$. From (42) and (36) we have

[TABLE]

and we deduce that there exists $k_{3}>0$ such that $0<k\leq k_{3}$ implies

[TABLE]

where $C_{i}$ is a constant depending only on $k_{3}$, $T$, and the derivatives of $u$ and $F$ up to order $3$ and $i$, respectively. From the inductions hypothesis (35) and inequality (4) we have

[TABLE]

and

[TABLE]

where $C$ is a constant depending only on $m$, $T$, and the derivatives of $u$ and $F$ up to order $r+2$ and $r$, respectively. Each $L^{n,q}_{i,\alpha_{i},\beta_{i}}$ being multilinear continuous, we deduce from (44)-(46) and the relation $\beta_{i}^{l}+\alpha_{i}^{l}\leq q-l+1$, for $l=1,\cdots,i$, that

[TABLE]

It follows by the triangle inequality that (43) for $q=2m+1$ yields

[TABLE]

for $n=m,m+1,\cdots,N-(m+1)-1$, where $C$ is a constant depending only on $p$, $T$, and the derivatives of $u$ and $F$ up to order $2m+4$ and $2m+2$, respectively . Passing to the norm in identity (41), we deduce from (38) and the last inequality that

[TABLE]

Otherwise, applying $D_{-}$ to (41), inequalities (44)-(46) and (47) yield

[TABLE]

for $n=m,m+1,\cdots,N-(m+1)-1$, where $C$ is a constant depending only on $p$, $T$, and the derivatives of $u$ and $F$ up to order $2m+5$ and $2m+3$, respectively. Therefore, passing to the norm in the identity obtained by applying $D_{-}$ to (41), we deduce from (41) and the last inequality that

[TABLE]

for $n=m,m+1,\cdots,N-(m+1)-1$, with the constant $C$ depending only on $p$, $T$, and the derivatives of $u$ and $F$ up to order $2m+5$ and $2m+3$, respectively. Inequalities (47) and (48) imply that the induction hypothesis is also true for $m+1$, and we deduce that (35) is true for each integer $m=0,1,...,p$.

Proof (Proof of Theorem 4.1)

We proceed by induction on $j=1,2,...,p+1$. The case $j=1$ is immediate from Lemma 1. Suppose that $\left\{u^{2j,n}\right\}_{n}^{N}$ approximates $u$ with order $2j$ of accuracy and satisfies (34), for an arbitrary $j$ such that $j\leq p$. We are going to prove that $\left\{u^{2j+2,n}\right\}_{n}^{N}$ approximates $u$ with order $2j+2$ of accuracy and (34) holds substituting $j$ by $j+1$.

From the induction hypothesis, $\left\{u^{2j,n}\right\}_{n}$ satisfies DCC. Because $\left\{u^{2j,n}\right\}_{n}$ and $\left\{\overline{u}^{2j,m}\right\}_{m}$ are computed from the same scheme DC2j, but for different time steps, $\left\{\overline{u}^{2j,m}\right\}_{m}$ also satisfies DCC. Therefore, as in 29, Theorem 3.1 applied to the approximation $\left\{u^{2j+2,n}\right\}_{n=0}^{j}$, built from $\left\{\overline{u}^{2j,m}\right\}_{m}$, yields

[TABLE]

where

[TABLE]

According to the DCC and the condition $u^{2j+2,0}=u(t_{0})=u_{0}$, we have

[TABLE]

By the triangle inequality and the DCC, the last two inequalities yield

[TABLE]

From the DCC on $\left\{u^{2j,n}\right\}_{n}$ and the inequality (49), Theorem 3.1 again implies that $\left\{u^{2j+2,n}\right\}_{n=0}^{N}$ approximates the exact solution $u$ with order $2j+2$ of accuracy. Therefore, it is enough to establish (34) for $j+1$, $j\leq p$. To this end we rewrite identity (24) as follows

[TABLE]

with

[TABLE]

where $\Theta^{2j+2,n}$ and $\sigma^{2j+2,n+1/2}$ are as in Theorem 3.1. Proceeding as in Lemma 1 and taking the finite difference formulae (5) and (6) into account, we can write

[TABLE]

where

[TABLE]

$C$ is a constant depending only on $p$, $T$, and the derivatives of $u$ and $F$. According to the inequality (34) from the induction hypothesis, we may write

[TABLE]

where

[TABLE]

Therefore, writing (50) as follows

[TABLE]

with

[TABLE]

the induction hypothesis and the reasoning from Lemma 1, substituting the functions $h$ and $g$, respectively, by $H$ and $G$, $\widehat{\Theta}^{2,n+1}$ by $\widehat{\Theta}^{2j+2,n+1}$, and $k^{2}$ by $k^{2j+2}$, yields

[TABLE]

for $m=0,1,...,p-j$ and $n=m+j-1,m+j,...,N-j-m$, where $C$ is a constant depending only on $p$, $T$, and the derivatives of $u$ and $F$ up to order $2(m+j+1)+1$ and $2(m+j)+1$, respectively. Inequality (34) holds for $\left\{u^{2j+2,n}\right\}_{n}$ by the triangle inequality from the last inequality.

We end this section by the following corollary that gives an important convergence property of the DC method. This property is useful for a time-stepping method to solve stiff and large dimensional differential equations arising from the space discretization of time-dependent PDEs.

Corollary 1

Suppose that the function $F$ is from $\mathbb{R}^{s}\rightarrow\mathbb{R}^{s}$, for a positive integer $s$, and satisfies the one-sided Lipschitz condition (22). Then, each approximate solution $\left\{u^{2j,n}\right\}_{n=0}^{N}$ from $DC2j$ satisfies the inequality

[TABLE]

where $C$ is a constant independent from any global Lipschitz constant on $F$, and either $k_{0}=2/\beta$ for $\beta>0$ or $k_{0}=+\infty$ for $\beta\leq 0$.

Proof

From the regularity assumption on $F$ and $u$ and the one sided-Lipschitz condition, we deduce from Theorem 4.1 that each $\left\{u^{2j,n}\right\}_{n=0}^{N}$, $j=1,2,\cdots$, satisfies DCC. Therefore, inequality (51) is immediate from the part (4) of Theorem 3.1. The constant $C$ depends only on the derivatives of $u$ up to order $2j+1$ and, according to (31)-(32) and the mean value theorem, on the bound of the Jacobian $F_{y}$ on the compact set $[0,T]\times\left\{y\in\mathbb{R}^{s}:|y|\leq R_{j}\right\}$.

Remark 6

The convergence property satisfied by the schemes $DC2j$ in Corollary 1 is in fact $B$-convergence (see, e.g., frank1981concept ; kraaijevanger1985b ) since the constant $C$ of the global error in (51) is independent from any global Lipschitz constant of the function $F$. Nevertheless, since in the definition of $B$-convergence the constant $C$ depends on high order derivatives of the exact solution $u$, the identity

[TABLE]

can make any requirement on the independence of the constant $C$ with respect to $F_{u}$ somewhat artificial. The numerical test on Bernoulli ODE in Section 6 gives an application of Corollary 1.

Remark 7

In practice, from part 4 of the proof of Theorem 3.1, the global error for an approximate solution of the IVP (1) under the one-sided Lipschitz condition (22) by a DC2j+2 method, $j=0,1,2,\cdots$, takes the form

[TABLE]

for $-2\leq\beta k<2$. The constant $C$ depends on the derivative of the exact solution $u$ of order $2j+3$ and can be very large in magnitude. However, if $\beta<0$ and $k$ is not too small, the factor $\left(\frac{2+\beta k}{2-\beta k}\right)^{n}$ is sufficiently small such that $C\left(\frac{2+\beta k}{2-\beta k}\right)^{n}<<1$, leading to very accurate approximate solutions for large time steps $k$. Nevertheless, independently of the sign of $\beta$, when $k$ is sufficiently small in the asymptotic region $k\mu<2$, where $\mu$ is the global Lipschitz constant of $F$, $\left(\frac{2+\beta k}{2-\beta k}\right)^{n}$ becomes closed to 1, for example when $n=j+1$, so that only $c_{2j+3}k^{2j+2}$ must dominate the constant $C$. Consequently, a non B-convergent method can be competitive with respect to a B-convergent one for sufficiently small time steps. This situation will be illustrated by the Bernoulli ODE in Section 6.

5 Absolute stability

In this section we prove the absolute stability of the DC schemes. The notion of absolute stability is introduced by Dahlquist MR0170477 to characterize methods able to solve stiff ODEs. Considering the following IVP,

[TABLE]

where $\lambda$ is a complex number, we have the following definition (see quarteroni2010 ; MR0170477 ):

Definition 2

A numerical method is said to be absolutely stable if the corresponding solution for the problem (53) for fixed $k>0$ and some $Re(\lambda)<0$ is such that

[TABLE]

The region of absolute stability of a numerical method is defined as the subset of the complex plane

[TABLE]

If $\mathcal{A}\cap\mathbb{C}_{-}=\mathbb{C}_{-}$, $\mathbb{C}_{-}=\left\{\lambda\in\mathbb{C}:Re(\lambda)<0\right\}$, the numerical method is said to be A-stable.

Before establishing absolute stability results for the deferred correction schemes (8) and (LABEL:a27)-(10), we recall the following result.

Lemma 2 ( see (tuenter2006frobenius, , formula (6)) )

Let $P_{m}$ be a polynomial of degree $m$ in one variable. Then the sum $\sum_{i=0}^{n}P_{m}(i)$ is a polynomial of degree $m+1$ in the variable $n$.

Lemma 3

Suppose that $F(t,u)=\lambda u$ and $u_{0}=1$ in the initial value problem (1), where $\lambda$ is a complex number with negative real part ($\lambda\in\mathbb{C}_{-}$). Then the corresponding approximate solutions from the schemes (8) and (LABEL:a27)-(10) can be written as follows

[TABLE]

where $P_{j}(n)$ is a polynomial of degree $j$ in the variable $n$.

Proof

We suppose that $\lambda k\neq-2$, otherwise we trivially have $u^{2j,n+1}=0$, for $n\geq j$. Since $F(t,u)=\lambda u$, we can rewrite (LABEL:a27) as follows

[TABLE]

where, according to formulae (2) and (3), we have

[TABLE]

and

[TABLE]

Combining the last three identities, we deduce that

[TABLE]

where $\alpha_{j,i}$ is affine in $\lambda k$. Under the hypothesis of the lemma, (8) matches the trapezoidal rule, and we have

[TABLE]

that is (56) is true for $j=0$. Suppose that (56) holds for an arbitrary integer $j\geq 0$. From (57) we have

[TABLE]

with $n\geq j+2$, and, substituting each $u^{2j+2,n+1+j-i}$ by the formula given by the induction hypothesis (56), we deduce that

[TABLE]

where

[TABLE]

It follows that

[TABLE]

It is clear that $Q_{j}(n)$ is a polynomial of degree $j$ in the variable $n$ as $P_{j}(n)$. Therefore, according to the Lemma 2, $\sum_{i=j+2}^{n}Q_{j}(i)$ is a polynomial of degree $(j+1)$ in the variable $n$. Whence,

[TABLE]

where

[TABLE]

is a polynomial of degree $j+1$ in the variable $n$. We then deduce by induction that the lemma is true for arbitrary non-negative integer $j$.

Theorem 5.1

Each of the deferred correction schemes (8) and (LABEL:a27)-(10) is A-stable.

Proof

From Lemma 3 we have, for $Re(\lambda k)<0$,

[TABLE]

since, under the condition $Re(\lambda k)<0$, we have $\left|\frac{2+\lambda k}{2-\lambda k}\right|<1$.

6 Numerical experiments

In this section we evaluate the accuracy and order of convergence of the schemes $DC2,DC4,\cdots,DC10$, implemented using the Scilab programming language. The starting values are computed using the scheme (11)-(12).

We choose six standard problems for the evaluation. The first problem concerns $B$-convergence by considering a Bernoulli equation. The second problem is about long term integration with an oscillatory solution of large amplitude. The four other problems are about stiffness. The third and fourth problems (B5 modified and E5, respectively) both involve complex eigenvalues of negative real parts, where the imaginary parts of the eigenvalues for the third problem have larger magnitudes while those from the fourth problem have smaller magnitudes. The fifth problem (Robertson) is nonlinear and stiff with real negative eigenvalues, and it also addresses B-convergence. The sixth problem is the van der Pol oscillator, which is stiff with arbitrary complex eigenvalues.

The first three problems have analytic solutions. For problems (61), (62) and (63) that do not have an analytic solution, we consider a small time step such that the approximate solutions with $DC6,\cdots,DC10$ are almost identical (to machine precision for problem (62)), and we choose one of the approximate solutions as reference solution.

For solutions $u=(u_{1},\cdots,u_{d})~{}:~{}[0,T]\rightarrow\mathbb{R}^{d}$, $1\leq d\leq 6$, the absolute error on the approximate solutions $\left\{u^{2j,n}\right\}_{0\leq n\leq N}$, $1\leq j\leq 5$, is computed with the norm

[TABLE]

For very large $N$ we extract solutions at $2\times 10^{6}$ or $3\times 10^{6}$ discrete times evenly spread over the interval $[0,T]$.

For a comparison of accuracy, we implement in Scilab the backward differentiation formulae (BDF) of order 2, 4 and 6, and the explicit Runge-Kutta (RK) of order 4. The implemented BDF are run with exact starting values for the first three problems that have analytic solutions, while for problems four and five the starting values are provided by the function stiff (implementing BDF with adaptive steps) of the solver ode from Scilab. For the van der Pol oscillator, the comparison of our DC methods is done only with the solutions from stiff and rkf from the solver ode. For each of the problems, except the first one, we give a table of absolute errors and orders of convergence for pairs of two consecutive time steps, for the approximate solutions with the DC methods. We denote by $k_{max}$ the maximal time step allowed to compute an approximate solution with the solver stiff or rkf (see enright1975comparing for a discussion on maximal time steps).

6.1 Bernoulli differential equation

[TABLE]

Table 2 gives the absolute error and the order of convergence for each pair of consecutive time steps, in the case of DC, BDF and RK4 methods. The dash for RK4 indicates that the method is unstable for the corresponding time steps.

This problem addresses $B$-convergence since the function $F$ is one-sided Lipschitz with $\beta=-0.1$, when positive solutions are considered. Moreover, the problem is strongly nonlinear with exponentially increasing magnitude of derivatives of the right side function $F$. Such derivatives of large magnitude generally limit the accuracy of high order methods that are not B-convergent. The one-sided Lipschitz constant being negative, in accordance with Corollary 1, DC methods provide very accurate approximate solutions for large time steps, and their accuracy increases with the order of the method. However, the convergence of the DC methods is suboptimal, due to the effect of the strong nonlinearity of the ODE. While $DC4$ and $DC6$ almost achieve their proper order for $k\leq 3\times 10^{-5}$, the order of convergence of $DC8$ and $DC10$ are not observed since these methods quickly achieve machine accuracy. In fact, DC10 achieves order 8.05 of convergence for $k=7.14\times 10^{-6}$ to $k=6.66\times 10^{-6}$. BDF methods are stable for large time steps, but they are less accurate than their corresponding DC methods. RK4 is completely unstable for $k\geq 2.03\times 10^{-3}$. For sufficiently small time steps in the asymptotic region, RK4 is more accurate than DC4 and any of the BDF methods, as stated in Remark 7, while DC6-10 achieve better accuracy.

6.2 Oscillatory problem hull1972comparing

[TABLE]

The exact solution is $u(t)=e^{\lambda\sin(t)}$. The original problem is set with $\lambda=1$ in hull1972comparing . The author in karouma2015class solved this problem with Runge-Kutta methods of orders 4 and 8, for $\lambda=2$ and $T=2580\pi$, to “illustrate the need of higher order methods when a long-term integration problem is considered”. Table 3 gives the absolute error and the order of convergence for each pair of consecutive time steps. The BDF methods are run only for the smallest time step. The solvers rkf and stiff use adaptive time stepping with a maximal time step $k_{max}=0.1$ and tolerances $rtol=100\times atol=10^{-10}$.

The magnitude of the exact solution $u(t)=e^{10\sin(t)}$ of the modified oscillatory problem is large, resulting in a relatively large absolute error obtained by the DC schemes (absolute errors of about $10^{-7}$ is possible for a good choice of stepsize). Moreover, the long term integration influences the accuracy of these schemes since they achieve absolute errors of about $10^{-9}$ when the solution interval is reduced to $[0,1000]$. Nevertheless, each DC scheme converges with its proper order. The DC methods are considerably more accurate than standard methods (both with fixed and variable stepsizes) which are inaccurate for this problem. For instance, for BDF2 and rkf, the solutions remain bounded with bounds close to the maximal amplitude of the exact solution but the phase of the oscillation is completely wrong.

6.3 Problem B5 modified enright1975comparing , stiff with complex eigenvalues of negative real parts and larger (in magnitude) imaginary parts

[TABLE]

This problem, originally set with $\alpha=100$, is an illustration of ODEs resulting from a semi-discretization by finite element methods of parabolic PDEs stewart1990avoiding . We choose $\alpha=5000$ to make the problem a little more difficult. Table 4 gives the absolute errors for the first component of the approximate solutions which is similar for the second component. The absolute errors for the others components quickly achieve machine precision. The solvers stiff and rkf are run with $k_{max}=2\times 10^{-5}$ and $atol=10\times rtol=10^{-15}$.

The imaginary parts of the Jacobian eigenvalues of the modified B5 problem are large. Even though the real parts of the eigenvalues are negative, we observe that smaller time steps are required by DC schemes to obtain accurate approximations. DC schemes achieve their proper order of convergence, but BDF methods perform better for this problem than DC schemes.

6.4 Problem E5enright1975comparing , stiff with complex eigenvalues of negative real parts and smaller (in magnitude) imaginary parts

[TABLE]

A reference solution is computed with $DC10$ for $k=10^{-3}$. The solution of this problem has small magnitude in $[1.618\times 10^{-3},1.76\times 10^{-3}]\times[0,1.46\times 10^{-10}]\times[0,8.27\times 10^{-12}]\times[0,1.38\times 10^{-10}]$ and the eigenvalues of the Jacobian matrix $dF(y)$ along the solution curve belong to the region $[-20490,3.68\times 10^{-12}]\times[-9.17\times 10^{-5},9.17\times 10^{-5}]$ of the complex plane. Table 5 gives the absolute errors and order of convergence for the four components of the approximate solutions. For BDF, RK4 and stiff, the absolute errors are provided only for the first component. The absolute error on the other components is smaller by 2 (RK4) to 5 (stiff) orders of magnitude, as we should expect from the magnitude of the solution components. The implemented BDF methods are run with starting values deduced from the solver stiff. The implemented RK4 is unstable for time steps $k\geq 2\times 10^{-4}$, and the absolute error is reported for $k=10^{-4}$ in table 5. The solver stiff is run with $k_{max}=10^{-3}$ and $rtol=10^{8}\times atol=10^{-15}$.

Imaginary parts of eigenvalues for the problem E5 are smaller, and larger time steps allow DC schemes to produce very accurate approximations, compared to the modified B5 problem. DC schemes perform better for this problem than BDF methods. They achieve their proper order of convergence but on a relatively small range of time steps, for higher order DC methods, since the solution is already very accurate for large time steps.

6.5 Robertson (1966) wanner1991solving , stiff with real negative eigenvalues

[TABLE]

This is one of the three problems considered as stiffest in wanner1991solving . We compute a reference solution with DC10 for the time step $k=1/6000$. The solution belongs to the region $[1.78\times 10^{-2},1.00]\times[0,3.58\times 10^{-5}]\times[0,0.983]$ and the eigenvalues of the Jacobian dF(y) along the solution curve belong to $[-9825.744,0]$. Table 6 gives absolute errors and orders of convergence of DC methods for each component of the solution. For other methods, we give only the maximal errors on the three components of the approximate solutions. The solver stiff is run with $k_{max}=1/600$ and $rtol=100\times atol=10^{-15}$. The solver rkf fails in solving this problem for various tolerances and $k_{max}$, and Scilab reported: “it is likely that rkf45 is inefficient for solving this problem”. The implemented BDF methods are run with starting values deduced from the solver stiff using the preceding tolerances.

The Robertson problem is stiff and addresses B-convergence since its Jacobian matrix has real negative eigenvalues with some having large magnitude. For this problem, DC schemes produce accurate approximate solutions even for large time steps, and high order DC methods can be avoided (DC6 is enough). The convergence is slow for $k>1/300$, but faster convergence happens for $k$ in the asymptotic region ($k<1/300$). The DC schemes perform better than BDF methods at equal order and time step. A comparison of the errors for $k=1/600$ suggests that the error constants might be 3 to 5 orders of magnitude smaller for DC than BDF methods.

6.6 van der Pol oscillator enright1975comparing ; shampine1981evaluation , stiff with arbitrary complex eigenvalues

[TABLE]

This problem was initially proposed for $T=1$ and $\mu=5$ in enright1975comparing . The actual version results from a suggestion by Shampine shampine1981evaluation . We compute a reference solution with DC8 for $k=1.875\times\times 10^{-6}$. The solution belong to the region $[-2,2.000073]\times[-1323.04,1231.35]$ of the real plane and the eigenvalues along the solution curve belong to the region $[-3000.29,1123.17]\times[-1158.48,1158.48]$ of the complex plane. Table 7 gives the absolute errors and orders of convergence. For rkf and stiff, we use $k_{max}=7.5\times 10^{-5}$ and $rtol=10$, $atol=10^{-16}$.

The van der Pol oscillator is stiff and the solution has a large magnitude. DC6 and DC8 reached their order of convergence. This shows that the DC strategy works well in spite of the fact that DC2 and DC4 would require much smaller time steps to produce reasonably accurate solutions. The order of convergence for DC10 is not observed, though the solutions obtained are accurate.

6.7 Discussion of the numerical results

In general, a careful assessment of the proof of Theorem 3.1 points out to the fact that, for a system with complex eigenvalues $\lambda=\lambda_{1}+i\lambda_{2}$, we only need a time step $k$ such that $k\,\max\left\{\lambda_{1},|\lambda_{2}|\right\}<2$ for a good accuracy (faster convergence happens when $-\lambda_{1}>>|\lambda_{2}|$). These situations are well illustrated by the test cases of Sections 6.3 and 6.4, where the required time step for accuracy is much smaller for modified B5 than E5. However, time steps $k$ such that $k\mu\simeq k|\lambda|<2$, $\mu\simeq\displaystyle\max_{0\leq t\leq T}\|d_{u}F\left(t,u(t)\right)\|$, is necessary for an asymptotic convergence with proper order. For example, in the case of the Bernoulli equation we have $\lambda\simeq-20000.1<0$ and $\mu=20000.1$. Large time steps provide accurate approximations (as expected from B-convergent methods), but asymptotic convergences are observed only for $k\mu<2$.

For the computational effort of the DC methods, we recall that to compute an approximate solution on discrete points $0=t_{0}<t_{1}<\cdots<t_{N}=T$, $DC2$ solves $N$ nonlinear systems while $DC2j$, $j\geq 2$, solves $j\times N$ systems. In the case of the Bernoulli equation, for example, $DC10$ achieves the maximal error of about $1.1\times 10^{-11}$ by solving approximately $5\times 10^{6}$ nonlinear systems while the maximal absolute error for $DC2$ is about $8.9\times 10^{-7}$ for $N=5\times 10^{6}$. We did not report any CPU time since our code is written in Scilab, an interpreted language. All methods that we implemented are consequently interpreted, while rkf and stiff provided with Scilab are compiled. Nevertheless, the main burden in implicit time-stepping solvers is the resolution of nonlinear systems, and we have shown that higher order DC methods give the most accurate approximations by solving fewer systems of equations. This gives a clue on the CPU time required and the efficiency of these methods. High order DC methods should be competitive in situations where using fully implicit methods is unavoidable.

7 Conclusions

We have presented a new approach of deferred correction methods for the numerical solution of general first order ordinary differential equations. Proofs for consistency, order of convergence and stability of the method are given, which rely on a recursive argument using a new deferred correction condition. The numerical experiments comply with the theory and show a high accuracy of the method and its satisfactory A-stable property and B-convergence. Globally, each DC scheme reaches its proper order of convergence and applies to any category of problem, providing accurate approximations for time steps not necessarily small. The accuracy of the DC schemes increases with the level of correction.