On the convergence of spectral deferred correction methods

TL;DR

This paper analyzes the convergence of spectral deferred correction (SDC) methods, revealing that the quadrature rule determines the order of accuracy, and establishes conditions for high-order convergence even with inconsistent solvers.

Contribution

It provides a new convergence analysis of SDC methods, emphasizing the role of quadrature rules and error sources, and extends the applicability to a broader class of solvers.

Findings

Quadrature rule determines the order of accuracy in SDC methods.

High-order accuracy can be achieved with inconsistent solvers under certain conditions.

Identifies three sources of errors in SDC methods: current error, previous error, and quadrature error.

Abstract

In this work we analyze the convergence properties of the Spectral Deferred Correction (SDC) method originally proposed by Dutt et al. [BIT, 40 (2000), pp. 241--266]. The framework for this high-order ordinary differential equation (ODE) solver is typically described wherein a low-order approximation (such as forward or backward Euler) is lifted to higher order accuracy by applying the same low-order method to an error equation and then adding in the resulting defect to correct the solution. Our focus is not on solving the error equation to increase the order of accuracy, but on rewriting the solver as an iterative Picard integral equation solver. In doing so, our chief finding is that it is not the low-order solver that picks up the order of accuracy with each correction, but it is the underlying quadrature rule of the right hand side function that is solely responsible for picking up additional orders of accuracy. Our proofs point to a total of three sources of errors that SDC methods carry: the error at the current time point, the error from the previous iterate, and the numerical integration error that comes from the total number of quadrature nodes used for integration. The second of these two sources of errors is what separates SDC methods from Picard integral equation methods; our findings indicate that as long as difference between the current and previous iterate always gets multiplied by at least a constant multiple of the time step size, then high-order accuracy can be found even if the underlying "solver" is inconsistent the underlying ODE. From this vantage, we solidify the prospects of extending spectral deferred correction methods to a larger class of solvers to which we present some examples.

Tables

Table 1. Table 1. Maximum size of the Lagrange polynomials max 1 ≤ n ≤ M ⁡ max ξ ∈ [ 0 , 1 ] ⁡ | ℓ n ​ ( ξ ) | subscript 1 𝑛 𝑀 subscript 𝜉 0 1 subscript ℓ 𝑛 𝜉 \max_{1\leq n\leq M}\max_{\xi\in[0,1]}|\ell_{n}(\xi)| for different quadrature points. The Gauss-Legendre, Gauss-Radau, and Gauss-Lobatto quadrature rules with M 𝑀 M points have degrees of precision 2 ​ M + 1 , 2 ​ M 2 𝑀 1 2 𝑀 2M+1,2M , and 2 ​ M − 1 2 𝑀 1 2M-1 , respectively.

	Type of Quadrature Points
M	Uniform	Chebyshev	Legendre	Gauss-Radau	Gauss-Lobatto
$2$	$1.000$	$1.207$	$1.366$	$1.500$	$1.000$
$3$	$1.000$	$1.244$	$1.479$	$1.558$	$1.000$
$4$	$1.056$	$1.257$	$1.527$	$1.578$	$1.000$
$5$	$1.152$	$1.263$	$1.551$	$1.586$	$1.000$
$6$	$1.257$	$1.266$	$1.566$	$1.591$	$1.000$
$7$	$1.362$	$1.268$	$1.575$	$1.594$	$1.000$
$8$	$1.663$	$1.269$	$1.581$	$1.596$	$1.000$
$9$	$2.550$	$1.270$	$1.585$	$1.597$	$1.000$
$10$	$4.028$	$1.271$	$1.588$	$1.598$	$1.000$
$11$	$6.506$	$1.271$	$1.590$	$1.599$	$1.000$
$12$	$10.963$	$1.271$	$1.592$	$1.599$	$1.000$
$13$	$18.340$	$1.272$	$1.594$	$1.600$	$1.000$
$14$	$32.060$	$1.272$	$1.595$	$1.600$	$1.000$
$15$	$54.998$	$1.272$	$1.596$	$1.600$	$1.000$
$16$	$98.531$	$1.272$	$1.596$	$1.600$	$1.000$
$17$	$172.176$	$1.272$	$1.597$	$1.601$	$1.000$
$18$	$313.675$	$1.272$	$1.597$	$1.601$	$1.000$
$19$	$556.491$	$1.273$	$1.598$	$1.601$	$1.000$
$20$	$1026.313$	$1.273$	$1.598$	$1.601$	$1.000$
$30$	$496210.554$	$1.273$	$1.600$	$1.602$	$1.000$
$50$	$208948162475.383$	$1.273$	$1.601$	$1.602$	$1.000$

Table 2. Table 2. Correction errors with M = 3 𝑀 3 M=3 uniformly spaced points.

Interval	Trapezoidal Rule		Forward Euler
	Error	$p=0$	Error	$p=0$
$[0,\frac{h}{2}]$	$\frac{h}{24}\left(e_1^{\left[p\right]}-2e_2^{\left[p\right]}+e_3^{\left[p\right]}\right)$	$𝒪\left(h^4\right)$	$\frac{h}{24}\left(7e_1^{\left[p\right]}-8e_2^{\left[p\right]}+e_3^{\left[p\right]}\right)$	$𝒪\left(h^3\right)$
$[\frac{h}{2},h]$	$\frac{h}{24}\left(e_1^{\left[p\right]}-2e_2^{\left[p\right]}+e_3^{\left[p\right]}\right)$	$𝒪\left(h^4\right)$	$\frac{h}{24}\left(e_1^{\left[p\right]}+4e_2^{\left[p\right]}-5e_3^{\left[p\right]}\right)$	$𝒪\left(h^3\right)$

Table 3. Table 3. Correction errors with M = 4 𝑀 4 M=4 uniformly spaced points.

Interval	Trapezoidal Rule		Forward Euler
	Error	$p=0$	Error	$p=0$
$[0,\frac{h}{3}]$	$\frac{h}{72}\left(3e_1^{\left[p\right]}-7e_2^{\left[p\right]}+5e_3^{\left[p\right]}-e_4^{\left[p\right]}\right)$	$𝒪\left(h^4\right)$	$\frac{h}{72}\left(15e_1^{\left[p\right]}-19e_2^{\left[p\right]}+5e_3^{\left[p\right]}-e_4^{\left[p\right]}\right)$	$𝒪\left(h^3\right)$
$[\frac{h}{3},\frac{2h}{3}]$	$\frac{h}{72}\left(e_1^{\left[p\right]}-e_2^{\left[p\right]}-e_3^{\left[p\right]}+e_4^{\left[p\right]}\right)$	$𝒪\left(h^4\right)$	$\frac{h}{72}\left(e_1^{\left[p\right]}+11e_2^{\left[p\right]}-13e_3^{\left[p\right]}+e_4^{\left[p\right]}\right)$	$𝒪\left(h^3\right)$
$[\frac{2h}{3},h]$	$\frac{h}{72}\left(e_1^{\left[p\right]}-5e_2^{\left[p\right]}+7e_3^{\left[p\right]}-3e_4^{\left[p\right]}\right)$	$𝒪\left(h^4\right)$	$\frac{h}{72}\left(e_1^{\left[p\right]}-5e_2^{\left[p\right]}-5e_3^{\left[p\right]}+9e_4^{\left[p\right]}\right)$	$𝒪\left(h^3\right)$

Table 4. Table 4. Correction errors with M = 4 𝑀 4 M=4 non-equispaced points. In this case, we find that both methods only pick up a single additional order of accuracy.

Interval	Trapezoidal Rule		Forward Euler
	Error	$p=0$	Error	$p=0$
$[0,\frac{h}{3}]$	$\frac{h}{162}\left(8e_1^{\left[p\right]}-27e_2^{\left[p\right]}+20e_3^{\left[p\right]}-e_4^{\left[p\right]}\right)$	$𝒪\left(h^3\right)$	$\frac{h}{162}\left(35e_1^{\left[p\right]}-54e_2^{\left[p\right]}+20e_3^{\left[p\right]}-e_4^{\left[p\right]}\right)$	$𝒪\left(h^3\right)$
$[\frac{h}{3},\frac{h}{2}]$	$\frac{h}{5184}\left(14e_1^{\left[p\right]}-27e_2^{\left[p\right]}+8e_3^{\left[p\right]}+5e_4^{\left[p\right]}\right)$	$𝒪\left(h^3\right)$	$\frac{h}{5184}\left(14e_1^{\left[p\right]}+405e_2^{\left[p\right]}-424e_3^{\left[p\right]}+5e_4^{\left[p\right]}\right)$	$𝒪\left(h^3\right)$
$[\frac{h}{2},h]$	$\frac{h}{192}\left(10e_1^{\left[p\right]}-81e_2^{\left[p\right]}+88e_3^{\left[p\right]}-17e_4^{\left[p\right]}\right)$	$𝒪\left(h^3\right)$	$\frac{h}{192}\left(10e_1^{\left[p\right]}-81e_2^{\left[p\right]}+40e_3^{\left[p\right]}+31e_4^{\left[p\right]}\right)$	$𝒪\left(h^3\right)$

Table 5. Table 5. Van der Pol oscillator. Here we present numerical results where we compare the implicit classical method defined in ( 24 ) as well as the modified semi-implicit SDC method defined in ( 66 ) against each other. Despite the fact that theory can show that the errors could be smaller for the modified method that relies solely on backward Euler time stepping embedded within Picard iteration, in this case the classical SISDC method based forward/backward Euler time stepping outperforms the other solver with its smaller error constants.

Mesh	SISDC	Order	Modified SISDC	Order
$4$	$2.24\times {10}^{-02}$	—	$6.45\times {10}^{-02}$	—
$8$	$6.06\times {10}^{-04}$	$5.21$	$2.84\times {10}^{-03}$	$4.51$
$16$	$4.11\times {10}^{-05}$	$3.88$	$1.91\times {10}^{-04}$	$3.89$
$32$	$3.44\times {10}^{-06}$	$3.58$	$1.46\times {10}^{-05}$	$3.71$
$64$	$2.56\times {10}^{-07}$	$3.75$	$1.01\times {10}^{-06}$	$3.85$
$128$	$1.78\times {10}^{-08}$	$3.85$	$6.68\times {10}^{-08}$	$3.92$
$256$	$1.17\times {10}^{-09}$	$3.93$	$4.29\times {10}^{-09}$	$3.96$
$512$	$7.26\times {10}^{-11}$	$4.01$	$2.69\times {10}^{-10}$	$3.99$

Equations

$$y^{\prime }=\frac{dy}{dt}$$

$$y(t)=y_0+{\int }_0^{t}f(y(s))ds,t>0.$$

$$|f(z)-f(w)|\le L|z-w|,$$

$$N=\{\begin{array}{ll}M-1& \text{if both endpoints are used},\\ M& \text{if only one endpoint is used},\\ M+1& \text{if neither endpoint is used}.\end{array}$$

$${\xi }_n^R=\{\begin{array}{ll}{\xi }_{n+1}& \text{if the left endpoint is included},\\ {\xi }_n& \text{if the left endpoint not included},\end{array}$$

$$w_{n,m}={\int }_{{\xi }_{n-1}^R}^{{\xi }_n^R}{\ell }_m(x)dx,n=1,2,\dots N,m=1,2,\dots M,$$

$${\ell }_m(x)=\frac{1}{c_m}\prod _{k=1,k\ne m}^M(x-{\xi }_k),c_m=\prod _{k=1,k\ne m}^M({\xi }_m-{\xi }_k).$$

$$\mathbf{Fully}\mathbf{Implicit}\mathbf{Collocation}:{\eta }_n={\eta }_{n-1}+h\sum _{m=1}^M{w}_{n,m}f({\eta }_m),n=1,2,\dots ,N,$$

$$y_n=y_{n-1}+{\int }_{t_{n-1}}^{t_n}f(y(t))dt,t_n={\xi }_n^{R}h,n=1,2,\dots N.$$

$$\mathbf{Picard}\mathbf{Iteration}:{\eta }_n^{[p+1]}={\eta }_{n-1}^{[p+1]}+h\sum _{m=1}^M{w}_{n,m}f({\eta }_m^{[p]}),n=1,2,\dots N.$$

$$Exp.\mathbf{SDC}:{\eta }_n^{[p+1]}={\eta }_{n-1}^{[p+1]}+h_n\left[f({\eta }_{n-1}^{[p+1]})-f({\eta }_{n-1}^{[p]})\right]+h\sum _{m=1}^M{w}_{n,m}f({\eta }_m^{[p]}),$$

$$Imp.\mathbf{SDC}:{\eta }_n^{[p+1]}={\eta }_{n-1}^{[p+1]}+h_n\left[f({\eta }_n^{[p+1]})-f({\eta }_n^{[p]})\right]+h\sum _{m=1}^M{w}_{n,m}f({\eta }_m^{[p]}).$$

$$I[f(\mathbf{\eta })](t):=\sum _{m=1}^{M}f({\eta }_m)l_m(t/h),f(\mathbf{\eta }):=\left(f({\eta }_1),f({\eta }_2),\cdots ,f({\eta }_M)\right).$$

$${\int }_{t_{n-1}}^{t_n}I[f(\mathbf{\eta })](t)dt=h\sum _{m=1}^M{w}_{n,m}f({\eta }_m)$$

$$\left|{\int }_{t_{n-1}}^{t_n}I[f(\mathbf{\eta })](t)-f(y(t))dt\right|\le h\Vert \mathbf{\eta }-\mathbf{y}\Vert W_{n}L+\frac{F}{M!}{h}^{M+1}.$$

$$\Vert \mathbf{e}\Vert :=\underset{1\le n\le M}{\max }\left|e_n\right|,\mathbf{e}=\left(e_1,e_2,\cdots ,e_M\right),$$

$$W_n:=\sum _{m=1}^M{\int }_{{\xi }_{n-1}^R}^{{\xi }_n^R}|l_m(\xi )|d\xi .$$

$$\left|{\int }_{t_{n-1}}^{t_n}I[f(\mathbf{\eta })](t)-f(y(t))dt\right|$$

$$+\left|{\int }_{t_{n-1}}^{t_n}I[f(\mathbf{y})](t)-f(y(t))dt\right|.$$

$$\left|{\int }_{t_{n-1}}^{t_n}I[f(\mathbf{\eta })](t)-I[f(\mathbf{y})](t)dt\right|$$

$$\le h\sum _{m=1}^M|{\omega }_{n,m}|\left|f({\eta }_m)-f(y_m)\right|$$

$$\le hL\sum _{m=1}^M|{\omega }_{n,m}|\left|{\eta }_m-y_m\right|$$

$$\le hL\Vert \mathbf{\eta }-\mathbf{y}\Vert \sum _{m=1}^M|{\omega }_{n,m}|.$$

$$\left|I[f(\mathbf{y})](t)-f(y(t))\right|=\left|\frac{(f\circ y{)}^{(M)}\left(z(t)\right)}{M!}\prod _{m=1}^M(t-t_m)dt\right|\le \frac{F}{M!}\prod _{m=1}^M|t-t_m|.$$

$$|a_n|\le A^n|a_0|+\frac{A^n-1}{A-1}B.$$

$$|a_n|\le A^n|a_0|+n{A}^{n-1}B$$

$$\frac{A^n-1}{A-1}=1+A+\cdots A^{n-1}\le n{A}^{n-1}$$

$$y^{\prime }=f(y),f(y)=f_I(y)+f_E(y),y(0)=y_0,$$

$${\eta }_n^{[0]}={\eta }_{n-1}^{[0]}+h_n{f}_I({\eta }_n^{[0]})+h_n{f}_E({\eta }_{n-1}^{[0]}),n=1,2,\dots N,$$

$$\mathbf{SISDC}:{\eta }_n^{[p+1]}$$

$$+h_n\left[f_E({\eta }_{n-1}^{[p+1]})-f_E({\eta }_{n-1}^{[p]})\right]+h\sum _{m=1}^M{w}_{n,m}f({\eta }_m^{[p]}),$$

$$|e_n^{[p+1]}|\le e^{Nh(2L_I+L_E)}|e_0|+C_{1}h\Vert {\mathbf{e}}^{[p]}\Vert +C_2{h}^{M+1},$$

$$C_1=2N{e}^{Nh(2L_I+L_E)}W,\text{and}{C}_2=2N{e}^{Nh(2L_I+L_E)}\frac{F}{M!}$$

Full text

1. Introduction

The spectral deferred correction (SDC) method defines a large class of ordinary differential equation (ODE) solvers that were originally introduced in 2000 by Dutt, Greengard and Rokhlin [12]. These type of methods are typically introduced by defining an error equation, and then repeatedly applying the same low-order solver to the error equation and adding the solution back into the current approximation in order to pick up an order of accuracy. This idea can be traced back to the work of Zadunaisky in 1976 [40], where the author sought out high-order solvers in order to reduce numerical roundoff errors for astronomical applications. Before introducing the classical SDC methods defined in [12], we stop here to point out some of the recent work that has been happening over the past two decades including [30, 28, 9, 26, 10, 7, 22]. We refer the interested reader to [32] for a nice list of references for the first of these last two decades. Here, we provide a sampling of some of the current topics of interest to the community.

Many variations of the original SDC method are being studied as part of an effort to expedite the convergence of the solver. The chief goal here is to reduce the total number of iterations required to obtain the same high-order accuracy of the original method. These methods include the option of using Krylov deferred correction methods [20, 21] as well as the multi-level SDC methods [27, 36]. The multi-level approach starts with a lower order interpolant and then successively increases the degree of the interpolant with each future sweep of the method. This has the primary advantage of decreasing the overall number of function evaluations that need to be conducted, but introduces additional complications involving the need to evaluate interpolating polynomials. In the same vein, higher order embedded integrators have been explored within the so-called integral deferred correction (IDC) framework [9, 10], where a moderate order solver (such as second or fourth-order Runge-Kutta method) is embedded inside a very high order SDC solver. With this framework, each successive correction increases the order by the same amount as that of the base solver. In addition, parallel in time solvers [5, 31, 8, 13] are being investigated as a mechanism to address the needs of modern high performance computing architectures, and adaptive time stepping options have been more recently investigated in [6]. This work is based upon the nice property that SDC methods naturally embed a lower order solver inside a higher order solver.

In addition to the above mentioned extensions, various semi-implicit formulations have, and are currently being explored. While the original solver was meant for classical non-linear ODEs, semi-implicit formulations have been conducted as early as 2003 [30] and are still an ongoing topic of research [29, 7]. The effect of the choice of correctors including second order semi-implicit solvers for the error equation has been conducted in [25], and an investigation into the efficiency of semi-implicit and multi-implicit spectral deferred correction methods for problems with varying temporal scales have been conducted in [26]. Related high-order operator splitting methods have been proposed in [15, 3, 4], where the focus is not on an implicit-explicit splitting, but rather on splitting the right hand side of the ODE into smaller systems that can be more readily inverted with each sweep of the solver.

Very recent work includes applications of the SDC framework to generate exponential integrators of arbitrary orders [2], exploring interesting LU decompositions of the implicit Butcher tableau on non-equispaced grids [38], further investigation into high-order operator splitting [11], a comparison of essentially non-oscillatory (ENO) versus piecewise parabolic method (PPM) coupled with SDC time integrators [23], and additional implicit-explicit (IMEX) splittings for fast-wave slow-wave splitting constructed from within the SDC framework [35].

It is not our aim to conduct a comprehensive review and comparison of all of these methods, rather it is our goal to present rigorous analysis of the original method that can be extended to these more complicated solvers. With that in mind, we now turn to a brief introduction of the spectral deferred correction framework, and in the process of doing so, we seek to directly compare this method with that of the Picard integral formulation of a numerical ODE solver.

1.1. Picard Iteration and the SDC Framework

We begin by giving a brief description of classical SDC methods. In doing so, we explain the differences between SDC and that of Picard iteration, which defines the cornerstone of the present work.

Classical SDC solvers are designed to solve initial value problems of the form

[TABLE]

where $y$ can be taken to be a vector of unkowns. The solution $y(t)$ can be expressed as an integral through formal integration:

[TABLE]

In this work, we assume that $f$ is Lipshitz continuous. That is, we assume

[TABLE]

for some constant $L\geq 0$ and all $z,w\in\mathbb{R}$. This is sufficient to guarantee existence and uniqueness for solutions of IVP (1), and produce rigorous numerical error bounds for SDC methods.

Consider a set of $M$ quadrature points $0\leq\xi_{1}<\ldots<\xi_{M}\leq 1$ that partition the unit interval into a total of $N$ disjoint subintervals, defined by

[TABLE]

We make this choice because a given quadrature rule may or may not include the endpoints of the interval, and this convention allows us to study Gaussian quadrature rules, uniformly spaced quadrature rules, Radau II quadrature rules and others all within the same context. With that in mind, we define the right endpoints $\xi^{R}_{n}$, for $n=0,1,\dots N-1$, of each the $N$ subintervals as

[TABLE]

and $\xi^{R}_{0}=0,\xi^{R}_{N}=1$ for the two boundary edge cases. Next, we define quadrature weights by

[TABLE]

where $\ell_{m}(x)$ is the Lagrange interpolating polynomial of degree at most $M-1$ corresponding to the quadrature point $\xi_{m}$:

[TABLE]

Once these weights are obtained, approximate integral solutions, say $\eta_{m}\approx y(\xi^{R}_{m}h)$ for $h>0$ and $m=0,1,\dots,N$, can be formed via

[TABLE]

whereas the exact solution $y_{m}:=y(\xi^{R}_{m}h)$ satisfies the exact integral

[TABLE]

By convention, $\eta_{0}:=y_{0}$ is known to high-order (because it comes from the previous time step), and $\eta_{N}\approx y(h)$ constitutes one “full” time step. Since each substep uses information from all substeps to construct the right-hand side, the solution is higher order, but also requires the solution of a nonlinear system of $M$ unknowns (one for each quadrature point) at each time step. Although this integrator has some very nice properties (e.g., it can be made to be symplectic and $L$-stable for suitably chosen quadrature points), it is not typically used in practice given the additional storage requirements and the larger matrices that need to be inverted for each time step. This is particularly relevant when it is used as the base solver for a partial differential equation, but even these bounds are being explored as a viable option for PDE solvers such as the discontinuous Galerkin method [33].

In place of the fully implicit collocation method, Picard iteration (with numerical quadrature) defines a solver by iterating on a current solution $\eta^{[p]}_{n}$, $p\in\mathbb{Z}_{\geq 0}$ and then creates a better approximation through

[TABLE]

Note that the current value $\eta^{[p+1]}_{0}:=\eta_{0}\approx y_{0}$ is a known value that is equal to the exact solution up to high-order. While this solver picks up a single order of accuracy with each correction, it has the unfortunate consequence of having a finite region of absolute stability.

The explicit spectral deferred correction framework is

[TABLE]

where $h_{n}=(\xi^{R}_{n}-\xi^{R}_{n-1})h$ is the length of the $n^{th}$ sub-interval. This solver also has a finite region of absolute stability.

Remark 1**.**

*Although traditional SDC methods were originally cast as a method that corrects a provisional solution by solving an error equation, some modern descriptions of the same solver identify Eqn. (9) as the base solver, which has the added benefit of pointing out a solid link between SDC methods and that of iterative Picard integral equation solvers. *

In order to construct methods that have more favorable regions of absolute stability for stiff problems, the implicit SDC framework exacts multiple backward Euler time steps through each iteration with

[TABLE]

Note that this framework allows for implicit and high-order solutions to be constructed with greater computational efficiency when compared to the fully implicit collocation solver defined in Eqn. (6) because smaller systems need to be inverted in order to take a single time step.

Remark 2**.**

It has been noted that the scaling in front of the $h_{n}$ term does not have impact on the order of accuracy [39]. It is our aim with this work to solidify that claim with rigorous numerical bounds, which we do for both the explicit and implicit solvers.

Before doing so, we point out an aside that is in common with all SDC solvers.

Remark 3**.**

If $\lim_{p\to\infty}\eta^{[p]}_{n}=\eta_{n}$ converges, then solutions to Eqns. (8), (9), and (10) converge to that of the fully implicit collocation method defined in Eqn. (6).

While some SDC methods work with a fixed number of iterates in order to obtain a desired order of accuracy, there are many examples in the literature where convergence of the SDC iterations to the fully implicit scheme is considered. For example, the work in [38] is wholly concerned with this convergence, and not the accuracy of the underlying method for fixed iterations. Moreover, for stiff problems, it is well understood that using a fixed number of iterations can lead to order reduction which negates the advantage of using SDC in the first place. In addition, the multilevel spectral deferred correction (MLSDC) methods also are typically iterated to a residual tolerance, since one cannot be sure that coarse level sweeps will provide enough increase in accuracy (or decrease in the residual) [36].

One key advantage of iterating an SDC method to convergence is that when this is done, the method inherits well known and desirable properties that the fully implicit collocation method enjoys. For example, Kuntzmann [24] and Butcher [1] separately point out that if a total of $M$ Gaussian quadrature points are used, then the fully implicit collocation method will have superconvergence order $\mathcal{O}(h^{2M})$. (For more details, we refer the interested reader to the excellent tomes of Hairer, Wanner et al [16, 17, 18]. For example, see Section II.7 of [16], Theorem 5.2 from [17], or Theorem 1.5 from [18].) In general, the maximum order of accuracy for the underlying solver with $M$ quadrature points is $\mathcal{O}(h^{2M})$ if they are Gauss-Legendre points, $\mathcal{O}(h^{2M-1})$ for the RadauIIA points, and $\mathcal{O}(h^{2M-2})$ for Gauss-Lobatto points. (The local truncation error is one order higher.) Uniform points have $\mathcal{O}(h^{M})$ order of convergence if $M$ is even, and $\mathcal{O}(h^{M+1})$ if $M$ is odd. The extra pickup in the order of accuracy is due to symmetry of the quadrature rule. (For example, $M=1$ points reproduces the so-called “midpoint” rule, $M=3$ reproduces Simpson’s rule, and $M=5$ yields Boole’s rule, each of which pick up an extra order of accuracy.)

1.2. An Outline of the Present Work

Despite the increasing popularity of spectral deferred correction solvers, very little work has been performed on convergence results for this large class of methods. The results that are currently in the literature [15, 9, 10, 19, 37] typically proceed via induction on the current order of the approximate solution, and they all hinge on solving the error equation, wherein the same low-order solver is applied and then a defect, or correction, is added back into the current solution in order to increase its overall order of accuracy. In other recent work [34], the authors consider SDC methods as fixed point iterations on a Neumann series expansion. There, the low-order method is viewed as an efficient preconditioner (in numerical linear algebra language), and the SDC iterations are thought of as simplified Newton iterations. Additionally, the work in [38] makes use of linear algebra techniques in order to optimize coefficients so that the method converges faster to the collocation solution for stiff problems.

In this work, we do not require the use of the error equation, nor do we work with any sort of defect such as that defined in [40], rather we instead focus on the Picard integral underpinnings inherent to all SDC methods. Our work solely uses fundamental numerical analysis tools: error estimates for numerical interpolation and integration. While these tools do rely on quadrature rules, our proofs are generic enough to accommodate any set of quadrature points, which are an ongoing discussion in terms of how to construct base solvers.

In this work, we prove rigorous error bounds for both implicit and explicit SDC methods, and in doing so, we expect the reader will find that these methods can be thought of as being built upon classical Picard iteration. Our results are applicable for general quadrature rules, but unlike the findings found in [37], where convergence is proven using the error equation, our work relies on the fundamental mechanics behind why the solver works. That is, we point out that the primary contributor to the order of accuracy of the solver lies within the integral of the residual, and not necessarily the application of any base solver to an error equation.

Indeed, our proofs follow in similar manner to that required to prove the Picard–Lindelöf theorem, but our proofs take into account numerical quadrature errors and do not rely on exact integration of the right hand side function $f(y)$. The primary differences between our proofs and that of the Picard-Lindelöf theorem are the following:

• •

Spectral deferred correction methods require the use of numerical quadrature to approximate the integrals presented in the Picard-Lindelöf theorem. Our error estimates take into account any errors resulting from quadrature rules.

• •

Each correction step in the implicit scheme defined in (10) requires a nonlinear inversion, whereas the Picard-Lindelöf theorem is typically proven using exact integration.

There are two main results in this work, one for explicit SDC methods, and one for implicit SDC methods. These are both found as corollaries to a single theorem on semi-implicit SDC. In each case, we produce rigorous error bounds that are applicable for generic quadrature rules. Furthermore, we find that there are a total of three sources of error that SDC methods carry: i) the error from the previous time step, ii) the error from the previous iterate, and iii) the error from the quadrature rule being used.

The outline of this paper is as follows. In Section 2, we present some necessary lemmas concerning error estimates for integrals of interpolants as well as some error estimates for sequences of inequalities that show up in our proofs. In Section 3 we present a convergence proof for the more general case of a semi-implicit SDC solver, and then immediately point out two corollaries that prove implicit and explicit SDC methods converge. In Section 4, we present results for an SDC method that makes use of a higher order base solver, the trapezoidal rule. In Section 5 we present some numerical results, where we compare explicit SDC methods with Picard iterative methods, we investigate modified implicit SDC methods, and we experiment with different semi-implicit formulations of SDC methods. Error estimates for all of these variants come from direct extensions of the proofs found in this work. Finally, some conclusions and suggestions for future work are drawn up in Section 6.

2. Preliminaries

We now point out a couple of important tools that we use to show that SDC solvers converge. Our aim is to focus on a single time step. Without loss of generality, from here on out we will focus on constructing a solution over the interval $[0,h]$, where $h$ is the time step size and we will assume that $\eta_{0}\approx y_{0}$ is a high-order approximation to the exact solution.

2.1. Error estimates for integrals of interpolants

If ${\boldsymbol{\eta}}=\left(\eta_{1},\eta_{2},\cdots,\eta_{M}\right)$ is a set of discrete values and $t\in[0,h]$ is a time interval we are interested in studying, we define the interpolation operator $I$ to be the projection onto the space of polynomials of degree at most $M-1$ via

[TABLE]

Note that this produces the integration identity

[TABLE]

after integrating (11) over a subinterval $[t_{n-1},t_{n}]:=[h\xi_{n-1}^{R},h\xi_{n}^{R}]$, and the weights are defined as in Eqn. (4).

Convergence results for both the explicit and the implicit SDC method (as well as Picard iteration) require the use of the following lemma.

Lemma 2.1**.**

Suppose that $f\circ y\in C^{M}([0,h])$, $\lVert\frac{d^{M}}{dt^{M}}\left(f\circ y\right)\rVert_{\infty}\leq F$, and that $f$ is Lipschitz continuous with Lipshitz constant $L$. Then we have the estimate

[TABLE]

where the discrete norm is defined by

[TABLE]

and the constant $W_{n}$ is defined by

[TABLE]

For a fixed quadrature rule, this constant is finite and independent of the function.

Proof.

Add and subtract the Lagrange interpolant $I[f({\boldsymbol{y}})](t)$ for $f\circ y$ inside the left hand side of Eqn. (13) and apply the triangle inequality:

[TABLE]

An estimate for the first of these two terms follows by linearity of the interpolation operator:

[TABLE]

The quadrature weights in this estimate are be bounded above by $|\omega_{n,m}|\leq\int_{\xi^{R}_{n-1}}^{\xi^{R}_{n}}\left|\ell_{m}(\xi)\right|d\xi$ and then summed over all $m$ to produce the constant $W_{n}$.

The second of the two integrals in (16) is a function solely of the smoothness of $f$ and the choice of the quadrature rule. That is, classical interpolation error estimates result in a bound on the $M^{th}$ derivative of $f\circ y$ through a single point $z(t)\in[0,h]$ that yields

[TABLE]

Because $|t-t_{m}|\leq h$ for each $m$, the result follows after integration. ∎

We stop to point out that due to the Runge-phenomenon, the coefficient $W_{n}$ defined in Equation (15) can become quite large if a large number of quadrature points are chosen for constructing the polynomial interpolants required for the SDC method. In Table 1, we demonstrate a few sample values when uniform, Chebyshev, Gauss-Legendre, Gauss-Radau, and Gauss-Lobatto quadrature nodes are used to construct the polynomial interpolants. Uniform quadrature points tend to start performing quite poorly in the teens; however even a small amount of points, say five or six, produces a high-order numerical method compared to other ODE solvers because in this regime the error constant is reasonable. The selection of quadrature points that minimizes this portion of the error constant is the Gauss-Lobatto nodes, but because convergence is found through refinement in $h$ rather than $p$, any of these points will produce a method that converges, provided the exact solution has a suitable degree of regularity.

2.2. Error estimates for sequences of inequalities

Finally, we require a second Lemma as well as a simple Corollary. Both of these are stated in [14] and their proofs are elementary.

Lemma 2.2**.**

If $\{a_{n}\}_{n\in\mathbb{Z}_{\geq 0}}$ is a sequence that satisfies $|a_{n}|\leq A|a_{n-1}|+B$ with $A\neq 1$, then

[TABLE]

Proof.

Recursively apply the inequality and sum the remaining finite geometric series. ∎

Corollary 2.3**.**

If $A>1$ and $\{a_{n}\}_{n\in\mathbb{Z}}$ is a sequence that satisfies $|a_{n}|\leq A|a_{n-1}|+B$, then

[TABLE]

for every $n$.

Proof.

By Lemma 2.2, the sequence satisfies Eqn. (19). We estimate the (finite) geometric series by

[TABLE]

because there are a total $n$ terms and each $A^{l}\leq A^{n-1}$ for $l=0,1,\dots,n-1$. ∎

With these preliminaries out of the way, we are now ready to state and prove our main result.

3. Convergence Results

In place of separately proving explicit and implicit results for (1), we instead consider an umbrella class of ODEs, defined through a semi-implicit formulation:

[TABLE]

where $f_{I}$ is to be treated implicitly, and $f_{E}$ is to be treated explicitly. We assume that both $f_{I}$ and $f_{E}$ have Lipshitz constants $L_{I}$ and $L_{E}$, respectively. In turn, this implies that $f$ has a Lipshitz constant of $L:=L_{I}+L_{E}$. In the case where $f_{I}\equiv 0$, we set $L_{I}=0$, and in the case where $f_{E}\equiv 0$, we set $L_{E}=0$.

The classical semi-implicit SDC (SISDC) method for (22) begins with a provisional solution, or initial guess $\eta^{[0]}_{n}\approx y(\xi_{n}h)$, that is typically defined with

[TABLE]

where $h_{n}=(\xi^{R}_{n}-\xi^{R}_{n-1})h$. This yield a first-order implicit-explicit (IMEX) predictor for the solution based upon a forward-backward Euler method. Our numerical (and analytical) results indicate the “predictor” step has little bearing on the overall order of accuracy of the solver. For example, it is possible to hold the solution constant for the initial iteration and still obtain high-order accuracy, albeit with one additional iteration.

The classical SISDC method [30] iterates on the provisional solution through

[TABLE]

where $\eta^{[p+1]}_{0}=\eta_{0}$ is a known quantity. This value is typically taken to be the result from the previous time step, and we assume that it is known to high-order accuracy. Our focus is on the local truncation error, in which case we assume that the error at time zero is non-zero. That is, we assume $e_{0}=\eta_{0}-y_{0}\neq 0$. Once the single step error is established, a global error can be directly found using textbook techniques. In the event where $f_{I}\equiv 0$, we end up with the explicit SDC method defined in (9), and when $f_{E}\equiv 0$, we end up the implicit SDC defined in (10).

We repeat that the collocation method defined in (6) requires simultaneously solving for each $\eta_{n}$ and is clearly more expensive than multiple applications of the backward Euler method found in (24), either on part of or the entire right hand side. In the event where $f_{I}\equiv 0$, then the method should be less expensive to run for a single time step, but the regions of absolute stability suffer [28, 29]. We also repeat that provided that if $\boldsymbol{\eta}^{[p]}$ converge as $p\to\infty$, then Eqn. (24) defines a solution to Eqn. (6). Proving which initial guesses converge to the fully implicit solver is beyond the scope of this work. Currently, our aim is to show that each correction step in the SDC framework picks up at least a single order of accuracy to the order predetermined by the quadrature rule.

3.1. Statement of the Main Result

Theorem 3.1**.**

The errors for a single step of the semi-implicit SDC method satisfy

[TABLE]

provided $hL_{I}<1/2$, and where $N$ is the number of intervals under consideration, $L:=L_{I}+L_{E}$ is the Lipschitz constant of $f$,

[TABLE]

are constants that depend only on $f$, the exact solution $y$, and the selection of quadrature points.

In subsection 3.3 we point out two corollaries to this result, one for implicit and one for explicit SDC, but before proving this theorem, we stop to point out an important observation that is applicable to any of the aforementioned methods.

Remark 4**.**

The statement of this theorem highlights that there are a total of three sources of error that SDC methods admit, which are ordered by appearance in the right hand side of (25):

1. (1)

the error at the current time step: $e_{0}=\eta_{0}-y_{0}$, 2. (2)

the error from the previous iterate (or predictor): ${\boldsymbol{e}}^{[p]}={\boldsymbol{\eta}}^{[p]}-{\bf y}$, and 3. (3)

the number of quadrature points, $M$.

The most important take-away is that because the error from the previous iterate, $\lVert{\bf e}^{[p]}\rVert$, gets multiplied by a factor of $h$, the error gets improved by one order of accuracy with each correction. Of course this order reaches a maximum order based upon the number of the quadrature points chosen, which can be seen in the third source of error. This can be improved by selecting quadrature points with superconvergence properties such as the Gaussian or Gauss-Lobatto quadrature points. Finally, please note that we make no comment about how the “previous” function values were found. This is intentionally done so, because we would like to focus our attention on the impact of what a single correction does the solution. In doing so, this permits the analysis to apply to parallel implementations of SDC methods where synchronizations between different correctors (threads) are seldom seen [5, 13].

3.2. Proof of the Main Result

Proof.

We subtract the exact equation (7) from (24) and find that the discrete error evolution equation is

[TABLE]

The last term in this summand can be estimated by appealing to Lemma 2.1 and observing

[TABLE]

We estimate the other terms by making use of their respective Lipshitz constants:

[TABLE]

The third line follows from the second by adding and subtracting $y_{n}$ to the inside each of the absolute values containing $|\eta_{n}^{[p+1]}-\eta_{n}^{[p]}|$ and $|\eta_{n-1}^{[p+1]}-\eta_{n-1}^{[p]}|$. Note that we also make use of the fact that $h_{n}\leq h$, although this too can be relaxed.

We continue by subtracting $hL_{I}|e_{n}^{[p+1]}|$ from both sides, dividing by $1-hL_{I}>0$, recognizing that $h_{n}<h$, and collecting the remaining terms involving the “explicit” portions:

[TABLE]

where we define $W:=\max_{1\leq n\leq N}\left(1+W_{n}L\right)$.

We make use of two separate estimates for $1/(1-hL_{I})$ to estimate the two terms found in the right hand side of (29). For the first term, we expand the geometric series and keep the first two terms:

[TABLE]

This is valid because $hL_{I}<1$. Additionally, $hL_{I}<1/2$, and therefore

[TABLE]

Together, these estimates imply that the first term can be estimated with

[TABLE]

For the second term, we have $1/(1-hL_{I})\leq 2$ for all $hL_{I}\in[0,1/2]$. This leads us to observe that

[TABLE]

Next, we appeal to Corollary 2.3 and make use of $A=e^{2hL_{I}}\left(1+hL_{E}\right)>1$ and $B=2\left(hW\lVert{\bf e}^{[p]}\rVert+\frac{F}{M!}h^{M+1}\right)$ to conclude that

[TABLE]

Since $1+hL_{E}\leq e^{hL_{E}}$ and $n\leq N$, we have the desired result. ∎

3.3. Corollaries of Main Result: Implicit and Explicit Error Estimates

With the general case proven in Theorem 3.1, we find results for both implicit, as well as explicit SDC solvers. An immediate corollary to Theorem 3.1 can be found by setting $f_{E}\equiv 0$, in which case $L_{I}$ becomes the Lipshitz constant for $f$, and the SISDC solver reduces to classical SDC with backward Euler defined in (10).

Corollary 3.2**.**

The errors for a single step of the implicit SDC method defined in (10) satisfy

[TABLE]

provided $h<1/(2L)$. The constants $C_{1}$ and $C_{2}$ depend only on the smoothness of $f$, the exact solution $y$, and the choice of quadrature points.

It is worth noting that the error estimate provided here is an asymptotic error estimate. That is, one key assumption that we have to make is that $h<1/(2L)$, which we do not have to make for the explicit case. Unfortunately, one key benefit of implicit solvers is that large time steps can be taken, in which case it is certainly possible that the solver does not obey this assumption. For these cases, a rigorous error estimate and analysis when $h>1/(2L)$ would make for an interesting result, which would be especially important for multiscale problems that contain large time scale separations. This observation is beyond the scope of the present work.

A related corollary for explicit solvers with tighter error bounds can be found. The result is the following.

Corollary 3.3**.**

The errors for a single step of the explicit SDC method defined in (9) satisfy

[TABLE]

where $N$ is the number of intervals under consideration, $L$ is the Lipschitz constant of $f$ for the ODE $y^{\prime}=f(y)$, and $C_{1}$ and $C_{2}$ are constants that depend only on $f$, the exact solution $y$, and the selection of quadrature points.

Proof.

Revisit the proof of Theorem 3.1 and replace the error estimate for $1/(1-hL_{I})\leq 2$ with 1 instead of 2. ∎

4. Convergence Proofs for Higher Order Base Solvers

We now consider the spectral deferred correction method with the implicit trapezoidal rule as its base solver:

[TABLE]

What makes this method interesting is that it picks up a total of two orders of accuracy with each correction. Note again, that in the absence of the terms that the factor $h_{n}/2$ multiplies, this method reduces to explicit Picard iteration, which picks up a single additional order of accuracy with each correction.

The result we focus on is the impact of each correction step, in which case any provisional solution may be used. In order to retain large regions of absolute stability reasonable methods include low-order implicit solvers such as backward Euler, or the second-order implicit trapezoidal (Crank-Nicholson) rule:

[TABLE]

In this section, we examine the interplay between the integral over the entire time interval, and the addition of extra integral terms that allows this solver to pick up additional orders of accuracy.

Let us define the exact value of the right hand side function as $f_{n}=f(y(t_{n}))$, the approximate value of the right hand side function as $f_{n}^{[p]}=f(\eta^{[p]}_{n})$, the local and global quadrature rules for integration over the subinterval $[t_{n-1},t_{n}]$ as

[TABLE]

These definitions allow us to compactly write the SDC method with the Trapezoidal rule defined in (37) to read

[TABLE]

Recall that the exact solution satisfies the integral equation (7), which we repeat:

[TABLE]

Theorem 4.1**.**

When coupled with the implicit Trapezoidal rule, the errors for a single step of the Spectral Deferred Correction method satisfy

[TABLE]

provided $hL<1$, $N$ is the number of intervals under consideration, $M$ is the number of points involved, $L$ is the Lipschitz constant of $f$, $F$ is an upper bound for the $M^{th}$ derivative of $f$, and the function $E^{[p]}(t)$ is the polynomial interpolant for the error in the the approximation of the right hand side function during the $p^{th}$ iterate defined by

[TABLE]

where $\Delta f^{[p]}_{m}:=f^{[p]}_{m}-f_{m}$ for each $m=1,2,\dots M$. The norm defined in (40) is the maximum absolute value of the second derivative of $E^{[p]}$:

[TABLE]

Proof.

We subtract the exact solution defined in Eqn. (7) from the SDC method based upon the trapezoidal rule defined in Eqn. (39) to end up with

[TABLE]

where $I_{n}:=H_{n}-\int_{t_{n-1}}^{t_{n}}f(y(t))dt$ is the difference between the high-order (discrete) integral and the exact integral of the right hand side.

We now estimate each of the three terms to the right of $e^{[p+1]}_{n-1}$ in Eqn. (43) separately, starting with the first term:

[TABLE]

which follows from using the Lipshitz continuity of $f$. The third term can be estimated by first recognizing that

[TABLE]

and then using Eqn. (18) (which requires assuming that $f\circ y\in C^{M}$) in order to yield

[TABLE]

where $F$ is any number that satisfies $\lVert\frac{d^{M}}{dt^{M}}\left(f\circ y\right)\rVert_{\infty}\leq F$.

Finally, we address the second, and most interesting term on the right hand side of Eqn. (43). The key observation comes from recognizing this term as the difference between a low-order (local) quadrature, $T_{n}$, and a high-order (global) quadrature $H_{n}$. Note that the exact integral of the polynomial interpolant $E^{[p]}$ over the subinterval $[t_{n-1},t_{n}]$ is

[TABLE]

and that

[TABLE]

is a low-order approximation to this integral. By textbook results, we have

[TABLE]

where $\xi_{n}$ is some number between $t_{n-1}$ and $t_{n}$. Together, this implies

[TABLE]

All together, inserting equations (44), (49), and (45) into (43), we have

[TABLE]

After replacing each $h_{n}\leq h$, rearranging, and assuming that $\frac{hL}{2}<1$, we have

[TABLE]

We now verify the two underscored inequalities involving the $1\pm hL/2$ terms asserted in (51). Identical to (32), we have

[TABLE]

after expanding the rational expression in terms of a geometric series, and assuming that $hL/2<1$ in order to retain convergence. Because $1+hL/2\leq 1+hL$, we have

[TABLE]

which verifies the first of the two underscored inequalities. For the second one, we need only assume that $hL<1$, which yields $1/(1-hL/2)<2$.

All together, we have

[TABLE]

which yields the desired result after appealing to Corollary 2.3 and using $n\leq N$.

∎

Remark 5**.**

The same decomposition for the error can be used for SDC coupled with Forward (or Backward) Euler. That is, identical to the decomposition found in (43), we can decompose the error as

[TABLE]

where $L_{n}$ and $L_{n}^{[p]}$ denote a “low-order” integral of the right hand side, but this time we use

[TABLE]

for Forward Euler, or instead

[TABLE]

for Backward Euler. The third term ${\bf III}$ is again $\mathcal{O}(h^{M+1})$, and the first term ${\bf I}$ can be bounded by a constant times $\left|e_{n}^{[p]}\right|+\left|e_{n-1}^{[p+1]}\right|$. The lack of additional order pickup can be found by observing that the second source of error instead satisfies

[TABLE]

where $\xi_{n}$ is some number between $t_{n-1}$ and $t_{n}$. In the following examples, we compare this term to that found from the Trapezoidal rule.

4.1. Examples

To illustrate the results of the Theorem presented in this section, we consider the linear test case

[TABLE]

and examine the errors produced by the SDC method when coupled with a higher order base solver. The order pickup for the SDC method can be found by examining the size of second source of error, defined in Equations (48) and (56), given by $h_{n}^{3}/12\left(E^{[p]}\right)^{\prime\prime}(\xi_{n})$ for the trapezoidal rule, and $h_{n}^{2}/2\left(E^{[p]}\right)^{\prime}(\xi_{n})$ for the Forward (or Backward) Euler method. Note that the form and size of this error is identical for either the Forward or Backward Euler base solver.

In the following examples, we work out the size of this term for a few different case studies. Taylor expansions are found by making use of the Maple software package.

4.1.1. $M=3$ uniformly spaced points

We first consider the results of using a total of $M=3$ equispaced quadrature points and look at the local Error over the subinterval $[t_{n-1},t_{n}]$ produced by the second term ${\bf II}$ defined in Equations (48) and (56). In the second set of columns in Table 2 look at the size of this term by writing out Taylor expansions for the first $p=0$ SDC iteration of a solver constructed by taking a Forward Euler provisional solution for the first time step. Note that in this case, each point satisfies $e^{[0]}_{n}=\mathcal{O}(h^{2})$ for each $n$, because the local truncation error (LTE) for Euler’s method is second order accurate. Therefore, the jump from $\mathcal{O}(h^{2})$ to $\mathcal{O}(h^{4})$ indicates that the Trapezoidal correction picks up an additional two orders of accuracy, whereas the jump from $\mathcal{O}(h^{2})$ to $\mathcal{O}(h^{3})$ picks up a single additional order of accuracy for the Euler base solver, which is consistent with the theory.

4.1.2. $M=4$ uniformly spaced points

We next consider the same problem with a total of $M=4$ equispaced quadrature points. Again, we look at the errors for the Trapezoidal method compared to the Forward Euler method after first constructing a provisional solution with the Forward Euler method. We again observe that the Trapezoidal rule improves the order of accuracy of the provisional solution by two factors, whereas Forward Euler (or likewise Backward Euler) only improves the order by one. Results for these quadrature points are presented in Table 3.

4.1.3. $M=4$ non-equispaced spaced points

Finally, we consider a case with a total of $M=4$ non-equispaced points. As an illustrative example, we consider the quadrature points $\xi_{1}=0,\xi_{2}=1/3,\xi_{3}=1/2$, and $\xi_{4}=1$. We find that the trapezoidal error only increases the order of the solver by one degree, which is consistent with the findings in [9], where the authors show that when the second-order Runge-Kutta method is used as a corrector, then the solver does not always pick up two orders of accuracy with each correction loop. Results for this problem are presented in Table 4.

5. Numerical results

The primary contribution of this work is to construct rigorous error estimates for classical SDC methods, and therefore we only include a couple of numerical results. An abundance of SDC examples applied to ordinary and partial differential equations can be found in the literature. One of our key goals here is to promulgate the fact that the primary source of high-order accuracy inherent in all SDC methods comes from its underlying Picard integral formulation, and not necessarily the “base solver,” and therefore we focus our results on nearby variations of classical SDC methods and demonstrate how classical SDC methods can be extended to produce related high-order solvers.

First, we introduce a comparison of errors (and stability regions) for explicit SDC vs. Picard iteration, second we explore modifications of the constant in front of an implicit SDC method, and finally we compare semi-implicit SDC and modified semi-implicit SDC solvers. For the sake of brevity the proposed modifications to SDC methods are not formally analyzed but straightforward extensions of the theorems presented in this work can be constructed to present formal error bounds for these methods. The numerical evidence presented here supports this claim.

5.1. A comparison of explicit SDC and Picard iteration

In this numerical example, we compare the errors and stability regions by applying the Picard iterative method defined in (8) to that of the explicit SDC defined in (9). In order to present an equal comparison of these two solvers, we consider identical initial guesses, or provisional solutions ${\boldsymbol{\eta}}^{[0]}$ based upon forward Euler time stepping and we work with uniform quadrature points for all of our test cases.

5.1.1. Errors for a linear test case

In Figure 1, we compare errors for the linear equation

[TABLE]

at a final time of $T=10$ with $\lambda=-2$ and $\lambda=-5$. Other orders and values of $\lambda$ show similar results where we observe slightly smaller error constants when using the Picard iterative method compared to the equivalent SDC method. This is consistent with the findings of Corollary 3.3, because the first error estimate

[TABLE]

could be tightened up to read

[TABLE]

which produces a smaller (provable) overall error for the Picard method when compared to the SDC method.

5.1.2. A comparison of regions of absolute stability for explicit methods

Next, we seek to compare regions of absolute stability for explicit SDC methods and their Picard iterative cousins. Here we observe that the stability regions are slightly improved when the “Euler term” in the time stepping is dropped from the SDC method. That is to say, we find that the Picard iterative methods generally have larger regions of absolute stability when compared to their SDC counterparts.

In order to demonstrate this, in Figure 2, we include a comparison of plots of the regions of absolute stability, defined by

[TABLE]

where $\rho(z)$ is the amplification factor for various quadrature rules for both of these methods, $z:=\lambda h$, and $\lambda$ is defined as in Eqn. (58). There, we compare methods of orders two through ten, all based on equispaced quadrature points, forward Euler time stepping for the provisional solution, and the minimum number of corrections required to reach the desired order of accuracy. (For example, the third order method uses two corrections and the fifth order method uses four corrections.) Similar to most explicit Runge-Kutta methods, we find that the regions of absolute stability increase as the order is increased, but there are also more function evaluations per time step.

5.2. Implicit SDC methods with modified backward Euler time steps

Here, we consider implicit SDC methods with a variable constant in front of the vanishing term:

[TABLE]

This same scaling has already been explored in [39] for SDC methods, but there the authors only consider the case where $1/2\leq\theta\leq 1$. With $\theta=0$, we have (explicit) Picard iteration (provided the provisional solution is modified), with $\theta=1$, we have the classical implicit SDC method, and with negative values of $\theta$ we have backward Euler solves on negative time steps; none of the these changes effect the overall order of accuracy, only the size of the error constant and the regions of absolute stability. Following the proof of the main Theorem in this work, we find the following result under a modified estimate for the time step size.

Theorem 5.1**.**

The errors for a single step of the the modified implicit SDC method defined in (60) satisfy

[TABLE]

*provided $|\theta|h<1/(2L)$. The constants *

[TABLE]

again depend only on the smoothness of $f$, the exact solution $y$, the choice of quadrature points, but this time they also depend on $\theta$.

Note that the value of $\theta=0$ minimizes the size of these constants, but it also produces poor regions of absolute stability. With $\theta\gg 1$ we have a method that is heavy handed on multiple backward Euler solves, and therefore it has a very large region of absolute stability, however these methods unfortunately introduce larger error constants. Small values of $\theta$ decrease these error constants, but they modify the regions of absolute stability to the point where they become finite and therefore undesirable as these are implicit methods.

5.2.1. A verification of high-order accuracy: The nonlinear pendulum problem

As a verification of the high-order accuracy of the solvers, we consider the equations of motion for a nonlinear pendulum:

[TABLE]

with appropriate initial conditions. If we perform the change of variables $y_{1}(t)=x(t)$ and $y_{2}(t)=x^{\prime}(t)$, we end up with the following first-order nonlinear system of equations that is equivalent to Eqn. (62):

[TABLE]

We consider initial conditions defined by $(y_{1}(0),y_{2}(0))=(0,1)$ and we integrate this problem to a final time of $T=10$. To compute a reference solution, we use MATLAB’s builtin ode45 with a relative tolerance of $10^{-12}$ and an absolute tolerance of $10^{-14}$.

We present convergence results for this problem in Figure 3. These results indicate that each method is indeed high-order independent of the value of $\theta$. For the sake of brevity, we only report results for methods with a total of $M=4$ equispaced quadrature points, but we also compare results for different number of corrections. (This underlying quadrature rule is also known as Simpson’s $3/8$ rule, which has a smaller error constant than Simpson’s rule that uses $M=3$ equispaced points, but it comes at the cost of an additional function evaluation.) We not only find that smaller values of $\theta$ produce smaller errors, which the theory supports, but we also demonstrate that more corrections for the methods with large $\theta$ values can help to decrease the errors.

5.2.2. A parameter study of regions of absolute stability for implicit methods

Given the results of the previous section, it would be tempting to want to set $\theta=0$ in order to reduce the total error. What is missing from this observation is an understanding of the regions of absolute stability. We now address this question. We find that small values of $\theta$ produce finite regions of absolute stability, and that large values of $\theta$ increase the regions of absolute stability (when compared to classical SDC methods) but they also increase the stiffness of each implicit solve. With that being said, larger time steps should be able to be taken, but as is pointed out in the previous section, larger errors are introduced. This reproduces the usual tradeoff between being able to take large time steps with large errors or being forced to take smaller time steps but at an increased computational cost.

In Figure 4 we present results for a third order method with various values of $\theta$. There we plot contour plots of the modulus of the amplification factor $|\rho(z)|$ in place of the boundary defined by $|\rho(z)|=1$ because if we were to plot the boundary then it would not be clear what parts are stable. In this sequence of images, we present results for various values of $\theta\in[-0.4,5]$. Larger values of $\theta$ such as $\theta=100$ look very similar to that of $\theta=5$. Tests on methods of other orders produce similar results involving transitions between finite and infinite regions of absolute stability as $\theta$ increases from [math]. The tradeoff between the size and shape of the stability regions for various quadrature rules is left for future work.

Before continuing, we stop to point out that diagonally implicit Runge-Kutta methods (on nonequispaced points) can be constructed from this very same framework. The point here is that because the term involving the difference $f(\eta^{[p+1]}_{n})-f(\eta^{[p]}_{n})$ in Eqn. (60) does not contribute to the overall order of accuracy; the scaling in front of this term can be modified so that each implicit solve uses the exact same time step, which would result in a singly diagonally implicit Runge-Kutta (SDIRK) method. This could be advantageous for easing the implementation of SDC methods in large scale code bases. In such a case the provisional solution would have to be modified in order to retain constant time steps for each stage in the solver. This would mean an extra correction or a (low-order) polynomial interpolation step would be necessary to not lose the starting accuracy found in the provisional solution, which would again modify the regions of absolute stability for solvers of various orders.

Similar modifications have recently been explored on nonequispaced points from a linear algebra perspective. In [38], the author makes use of this vantage and optimizes their solvers by modifying the coefficients in the fixed point iteration matrices. It is pointed out that there are a number of items that could be optimized, such as the spectral radius of the solver (in order to optimize the convergence rate of the sweeps), the matrix norm (for the purposes of reducing the error under the assumption of a small number of sweeps), the error at the final time, the average reduction factor in each sweep block (for the purposes of adaptively choosing the number of sweeps, which could include flexible or greedy sweeps), and so on. Even though SDC methods are a subset of Runge-Kutta methods, and all of these options can be found by looking at this more general class of methods, one key advantage SDC methods enjoy is they do not typically sacrifice the difficult order conditions that a more generic RK method would have to address. At the same time, SDC methods still have ample levers to tune for optimization purposes.

5.3. Modified semi-implicit SDC methods

Recall that the semi-implicit SDC (SISDC) method begins with a partition of the right hand side into two functions $f_{I}$ and $f_{E}$ via

[TABLE]

and then they apply a forward Euler/backward Euler (FE/BE) pair to the right hand side defined in (24):

[TABLE]

With a straightforward extension the results from the present work point out that high-order accuracy can be achieved where there are no forward Euler time steps on the explicit term $f_{E}(y)$. That is, we propose examining the following modified SISDC method:

[TABLE]

Note that the “base solver” for this method is not even consistent with the underlying ODE. This serves as another example of how the findings from this work permit modifications to classical SDC methods in order to produce nearby variations. As an additional benefit, this reduces the computational coding complexity by asking the user to only define $f$ and $f_{I}$ as opposed to $f$, $f_{I}$ and $f_{E}$. This type of modification can be readily found after understanding the source of high-order accuracy inherent to the SDC framework.

Given that the iterative Picard methods demonstrate smaller errors by dropping the forward Euler terms in the right hand of the iterations from the SDC solver, one might expect that this method has better accuracy than its SISDC parent. In the event when $f_{I}=0$ (or is small), this would be true because in that case we would be comparing explicit SDC to Picard iteration, and we have already shown that those errors are smaller for some problems. However, we will shortly see that this is not necessarily the case. Even though this modification does not affect the overall order of the solver, we will show that for the following test case it does not improve the total overall error of the solver. With that being said, we believe it is still important to understand the source of the overall order of the SDC solvers, because only then can new methods be developed from the existing framework.

5.3.1. Van der Pol’s equation

As a prototypical IMEX example, we include results for Van der Pol’s equation:

[TABLE]

with appropriate initial conditions. After making the usual transformation of $y_{1}(t)=x(t)$, $y_{2}(t)=\mu x^{\prime}(t)$, and rescaling time through $t\rightarrow t/\mu$, we have the following system of differential equations [30, 25]:

[TABLE]

In an IMEX setting, this problem is typically split into $f_{E}(y)=\left(y_{2},0\right)$ and $f_{I}=\left(0,\left(-y_{1}+(1-y_{1}^{2})y_{2}\right)/\epsilon\right)$ as an effort to account for the stiffness as $\epsilon\to 0$. For this problem we only seek to verify the high-order accuracy of the classical semi-implicit SDC method defined in (24), denoted by SISDC, as well as the modified solver defined in (66), denoted by “modified SISDC.” With this aim in mind, we set $\epsilon=1$ so that the equations remain non-stiff, and we integrate to a long final time of $T=4$. The initial conditions are the same as those found in an example in [25], which are $y_{1}(0)=2$ and $y_{2}(0)=-0.666666654321$. In Table 5, we compare a convergence study for the fourth-order versions of these two methods where we use a total of four equispaced quadrature points, a provisional solution defined by the split forward/backward Euler method, as well as three corrections in the solver. For this problem, we find that the classical SISDC method has slightly smaller errors, despite what theory might otherwise predict we could observe. For problems where $f_{I}$ is negligible or small, the modified method should outperform the SISDC solver. Other values for $\epsilon$ produce similar findings, where the usual order reduction can be found as $\epsilon$ approaches zero. In other cases, the two solvers have similar behavior, and both show high order accuracy for large values of $\epsilon$. For brevity, these other results are omitted.

6. Conclusions

In this work we present rigorous error bounds for both explicit and implicit spectral deferred correction methods. Unlike most presentations that introduce SDC methods as a method that iteratively corrects provisional solutions by solving an error equation, our work hinges on the fact that the basic solver can be recast as a variation on Picard iteration. This observation allows new SDC methods to be developed through modifications of the (forward or backward) Euler part of the iterative procedure. In addition, we present some analysis for SDC methods constructed with higher order base solvers. In the numerical results section we present some sample variations that serve to indicate that the choice of the base solver need not be consistent with the underlying ODE in order to obtain a method that converges. That is, our findings indicate that it is not important to use the same low-order solver for each correction step because the desired high-order accuracy can be found in the integral of the residual. However, the choice of the base solver certainly has an impact on the overall scheme. For example, up to the degree of precision of the underlying quadrature rule, the choice of the forward or backward Euler method or even an inconsistent base solver leads to a single pickup on the order of accuracy of the solver with each correction step, whereas the choice of a trapezoidal rule for a base solver yields two orders of pickup with each correction step. For stiff problems, an implicit method constructed with the backward Euler method (or some variation of it) is certainly preferable due to the larger regions of absolute stability. Future work involves further analysis of embedded high-order base solvers as well as exploring further modifications of the solver to modify regions of absolute stability of existing solvers for explicit, implicit, and semi-implicit SDC methods.

Acknowledgements. We would like to thank the anonymous referees for their thoughtful comments, suggestions that improved the quality of this manuscript, and recommendations for future research. The work of D.C. Seal was supported by the Naval Academy Research Council.