Differential Equations
Mathwrist White Paper Series

(October 5, 2023)

Abstract

This document presents a technical overview of solving differential equations in Mathwrist’s C++ Numerical Programming Library (NPL). The problem class we are targetting includes initial value problems for ordinary differential equations (ODE) and initial value/boundary problems for 1-d parabolic linear partial differential equations (PDE).

1 Ordinary Differential Equation

1.1 Introduction

For a smooth function $y(t)$ governed by an ordinary differential equation (ODE)

y^{\prime}(t)=f(t,y(t)),t\in[0,T]

, where $f(t,y(t))$ is continuous and satisfies Lipschitz condition wrt $y$ , then $y(t)$ has a unique solution $y(t)$ for $0\leq t\leq T$ given initial value $y(0)=y_{0}$ . See [1] for the details of Lipschitz condition, proofs and various numerical methods.

To numerically solve an ODE initial value problem, we generate the approximation $u_{i}$ of $y(t_{i})$ over a sequence of time steps $t_{i},i=0,\cdots,N$ with $u_{0}=y_{0}$ . Explicit one step methods have the general form,

u_{i+1}=u_{i}+h\phi(\cdot,u_{i}),h=t_{i+1}-t_{i}

, where the evaluation of $\phi(\cdot,u_{i})$ depends only on the approximation value $u_{i}$ at the previous step $t_{i}$ . The simplest one-step method is Euler’s method, which has global relative error $\mathcal{O}(h)$ with step size $h$ .

The family of Runge-Kutta methods are high order accurate 1-step methods, which require multiple evaluations of $f(t_{i}+\Delta t,u_{i}+k)$ for some $\Delta t\in[0,h]$ and some bump amount $k$ . Evaluation at a mid-way point $t_{i}+\Delta t$ could be infeasible in practice. For example, $f(t,y(t))$ is an implicit function and can be computed only on discretized grid points $[t_{i},t_{i+1},\cdots]$ . If $f(t_{i}+\Delta t,\cdot)$ indeed can be easily computed, Runge-Kutta methods usually perform better than multi-step methods.

Multi-step methods can also produce high order accuracy. A $m$ -step multi-step method has the general form

	$\displaystyle u_{i+1}$	$\displaystyle=$	$\displaystyle\alpha_{m-1}u_{i}+\alpha_{m-2}u_{i-1}+\cdots+\alpha_{0}u_{i+1-m}+$
			$\displaystyle h\left[\beta_{m}f(t_{i+1},u_{i+1})+\beta_{m-1}f(t_{i},u_{i})+% \cdots+\beta_{0}f(t_{i+1-m},u_{i+1-m})\right]$

If coefficient $\beta_{m}=0$ , the method is explicit or otherwise implicit. Explicit and implicit approaches can be combined to have the flavor of prediction and correction, i.e. $u_{i+1}$ is first computed from the explicit method and then substituted to the implicit formula. Prediction and correction can be repeated until some error control tolerance is satisfied. Alternatively, prediction and correction can be used as error control to construct variable step-size methods.

The family of Adams-Bashforth methods and Adams-Moulton methods are all multi-step methods. The advantage of multi-step methods is that $f(t,y(t))$ is computed exactly on grid points and is computed only once at each $t_{i}$ . However, multi-step methods often face stability issue when $h$ is not fine enough. For smooth problems though, multi-step methods could be excellent choices.

In the same spirit of using extrapolation to increase the accuracy of numerical integration, extrapolation can be used on top of the Euler’s method to increase the accuracy. This is a very powerful technique that is easy to implement and has the advantage of error control. See [1], section 5.8 for details.

In Mathwrist NPL, we choose an appropriate explicit or implicit ODE method as the base method and then increase the accuracy order of the base method within an Iterative Deferred Correction (IDeC) framework. Taking Euler’s method as the base method for example, our IDeC ODE solvers overall require less number evaluation of $y^{\prime}(t)$ than that of the extrapolation method to achieve the same level of accuracy.

1.2 Iterative Deferred Correction (IDeC)

The idea of deferred correction on a low order ODE method to obtain high order accuracy can be traced back to [2] and [3]. Convergence proof and error order analysis are discussed in [4]. We shall call the low order ODE method the “base method”. IDeC gives one extra order accuracy on top of the base method almost for free. Applying multiple error correction iterations could then raise the accuracy order to a desired level.

1.2.1 Classic IDeC Method

Consider the $i$ -th step of Euler’s method $u^{(1)}_{i}=u^{(1)}_{i-1}+hf(t_{i-1},u^{(1)}_{i-1})$ to solve the initial value problem

y^{\prime}(t)=f(t,y(t)),0\leq t\leq T,y(0)=y_{0}

(1)

Here we use the notation $u^{(1)}$ to indicate that it is a first order approximation of $y(t)$ . Now construct a smooth interpolation function $L(t)$ passing through $p+1$ number of points $\{u^{(1)}_{i-p},\cdots,u^{(1)}_{i}\}$ that have been computed by Euler’s method. The error function $e(t)=y(t)-L(t)$ then satisfies the ODE

e^{\prime}(t)=f(t,y(t))-L^{\prime}(t)

(2)

Writing the ODE in terms of $e(t)$ and $L(t)$ using the relation $y(t)=L(t)+e(t)$ , the ODE is $e^{\prime}(t)=g(t,e(t))$ , where $g(t,e)=f(t,L(t)+e)-L^{\prime}(t)$ . Given initial value $e(0)=0$ , we can again solve $e(t)$ using Euler’s method

e_{i+1}=e_{i}+hg(t_{i},e_{i})=e_{i}+h\left[f(t_{i},u^{(1)}_{i}+e_{i})-L^{% \prime}(t_{i})\right]

(3)

Since Euler’s method produces order $\mathcal{O}(h)$ relative error, Euler approximation on $e(t)$ then has error order $\mathcal{O}(h^{2})$ . So adding correction $e_{i}$ to $u^{(1)}_{i}$ produces a second order approximation,

u^{(2)}_{i}=u^{(1)}_{i}+e_{i}

This process can be repeated to raise the accuracy order. In Mathwrist NPL, we use Lagrange interpolation polynomial to construct $L(t)$ . The degree $p$ of the interpolation polynomial is restricted to be $3<p<7$ for practical reasons.

Let $q$ be the number of error correction iterations. We restrict $q<p$ . Extra error correction iterations generally do not increase the accuracy but only adding numerical noise. Our explicit IDeC solvers also allow user to choose the base method to be either the explicit Euler’s method or the Adams-Bashforth 3-step method.

1.2.2 Spectral IDeC Method

We can rewrite ODE initial value problem (1) in an integral equation form,

y(t)=y_{0}+\int_{0}^{t}f(s,y(s))ds

(4)

Assume we have obtained an order- $n$ approximation $u^{(n)}(t)$ of $y(t)$ . An accuracy measure of this approximation can be defined by the residual function,

\epsilon(t)=y_{0}+\int_{0}^{t}f(s,u^{(n)}(s))ds-u^{(n)}(t)

(5)

Substituting the true error function $e(t)=y(t)-u^{(n)}(t)$ back to equation (4), we have

u^{(n)}(t)+e(t)=y_{0}+\int_{0}^{t}f(s,u^{(n)}(s)+e(s))ds

Then using the definition (5) to obtain a Picard integral equation

e(t)=\int_{0}^{t}\left[f(s,u^{(n)}(s)+e(s))-f(s,u^{(n)}(s))\right]ds+\epsilon(t)

(6)

Let $g(t,e)=f(t,u^{(n)}(t)+e(t))-f(t,u^{(n)}(t))$ and simplify equation (6) to

e(t)=\int_{0}^{t}g(s,e(s))ds+\epsilon(t)

(7)

Now for knot points $t_{i}\in[0,h],i=0,\cdots,m,t_{0}=0\text{ and }t_{m}=h$ , we derive the explicit Euler scheme for the Picard integral equation (7)

e_{i+1}=e_{i}+g(t_{i},e_{i})\Delta t_{i}+\epsilon_{i+1}-\epsilon_{i}

(8)

The spectral IDeC method uses Gaussian quadrature to compute $\epsilon_{i}$ from equation (5). In Mathwrist NPL, we use a Legendre polynomial $L(x(t))$ of degree $p$ to compute the quadrature, where function $x(t)$ mapps $t\in[0,h]$ to $x\in[-1,1]$ . This choice has $2p-2$ order of accuracy when computing $\epsilon(t)$ .

In the spectral IDeC case, we restrict $3<p<12$ and $q<p$ . Same as classic IDeC, the base method could be either the explicit Euler’s method or the Adams-Bashforth 3-step method.

1.2.3 Implicit IDeC

Explicit ODE methods need to use small enough step size $h$ in order for the numerical solution to be stable. Implicit methods could be used in IDeC to address the stability issue. In fact, for a well-posed linear ODE problem, implicit methods could be unconditionally stable.

A n-dimensional system of ODE is linear if it is in a general form,

y^{\prime}(t)=\mathbf{A}(t)y(t)+s(t)

(9)

, where $\mathbf{A}(t)$ is a $n\times n$ matrix that only depends on $t$ . $s(t)$ is a n-dimensional external source function. When $s(t)=0$ , the system is homogeneous. System (9) is stable when all eigen values of $\mathbf{A}(t)$ are negative. We shall assume the problems to be solved are well posed in stability.

Let $\mathbf{y}$ and $\mathbf{s}$ denote vector values of $y(t)$ and $s(t)$ . Using implicit Euler as the base method of a classic IDeC, we have

\displaystyle\mathbf{y}_{i+1}=\mathbf{y}_{i}+h\mathbf{y}^{\prime}(t_{i+1})=% \mathbf{y}_{i}+h\left[\mathbf{A}_{i+1}\mathbf{y}_{i+1}+\mathbf{s}_{i+1}\right]

, which is to solve the linear system

\displaystyle\left[\mathbf{I}-h\mathbf{A}_{i+1}\right]\mathbf{y}_{i+1}=\mathbf% {y}_{i}+h\mathbf{s}_{i+1}

It can be shown that if all eigen values of $\mathbf{A}_{i+1}$ are negative, the implicit Euler method is unconditionally stable. In Mathwrist NPL, we offer the combination of classic IDeC or spectral IDeC with implicit Euler as the base method. This is particularly useful to solve stiff ODE problems.

2 One Dimensional Parabolic PDE

2.1 Introduction

We consider a two dimensional smooth function $v(t,x)$ in a rectangle area $x\in(a,b),t\in(0,T)$ that is governed by a general form of 1-d linear parabolic PDE

\frac{\partial v}{\partial t}=f(t,x)\frac{\partial^{2}v}{\partial x^{2}}+g(t,x% )\frac{\partial v}{\partial x}+h(t,x)v+s(t,x)

(10)

. Functions $f(t,x)$ , $g(t,x)$ and $h(t,x)$ at the right hand side are coefficient functions. For well-posed problems, $f(t,x)>0$ . Function $s(t,x)$ at the end of the equation is the external source function. A PDE is called homogeneous when $s(t,x)=0$ , or inhomogeneous otherwise.

We often refer $t$ as the temporal variable and $x$ as the spatial or state variable. 1-d parabolic PDE refers to the fact there is only one state variable $x$ . In a linear parabolic PDE, all coefficient functions and the source function only depend on $t$ and $x$ . We will also use partial derivative notations $v_{t},v_{x},v_{xx}$ interchangeably with notations in (10) for convenience.

We could introduce a partial derivative operator $\mathcal{A}(v)$ wrt spatial variable and write the PDE as an ODE form

	$\displaystyle v_{t}$	$\displaystyle=$	$\displaystyle\mathcal{A}(v)+s(t,x)$		(11)
	$\displaystyle\mathcal{A}(v)$	$\displaystyle=$	$\displaystyle f(t,x)v_{xx}+g(t,x)v_{x}+h(t,x)v$		(12)

For the numerical methods that we will discuss in the following sections, operator $\mathcal{A}(v)$ translates to a matrix-vector multiplication $\mathbf{A}\mathbf{v}$ . For example, $\mathbf{A}$ is constructed as the finite difference coefficient matrix by a finite difference scheme and $\mathbf{v}$ is the vector of $v(t,x)$ approximated at discretized node points.

For parabolic PDEs, the appropriate boundary conditions [8] are listed in table 1. Effectively, the Robin condition is a mixture of Dirichlet and Neumann conditions. Other boundary conditions lead to no solution.

Given an initial value condition $v(0,x)=u(x)$ and appropriate boundary conditions, a parabolic PDE has a unique and stable solution. The solution of $v(t,x)$ is immediately smoothed out for any $t>0$ even though the initial value $v(0,x)$ could be non-smooth.

Table 1: Appropriate Boundary Conditions for Parabolic PDE

Dirichlet	$\begin{aligned} v(t,a)&=&u_{1}(t)\cr v(t,b)&=&u_{2}(t)\end{aligned}$
Neumann	$\begin{aligned} v_{x}(t,a)&=&u_{1}(t)\cr v_{x}(t,b)&=&u_{2}(t)\end{aligned}$
Robin	$\begin{aligned} \alpha_{1}v(t,a)+\beta_{1}v_{x}(t,a)&=&u_{1}(t)\cr\alpha_{2}v(% t,b)+\beta_{2}v_{x}(t,b)&=&u_{2}(t)\end{aligned}$

2.2 Crank-Nicolson Method

Rearrange the PDE form (10) to

v_{t}=w(t,x,v,v_{x},v_{xx})+s(t,x)

(13)

, where we express the 0-th to 2nd order spatial derivative terms in a new function $w(t,x,v,v_{x},v_{xx})$ . Consider a mesh $\{t_{l=0,\cdots,N}\}\times\{x_{j=0,\cdots,J}\}$ that is uniformly discretized in the spatial dimension, e.g. $x_{j}=a+i\frac{b-a}{J}$ . Denote the function values at grid point $(t_{l},x_{j})$ as $f_{l,j}:=f(t_{l},x_{j})$ , $g_{l,j}:=g(t_{l},x_{j})$ , $h_{l,j}:=h(t_{l},x_{j})$ , $s_{l,j}:=s(t_{l},x_{j})$ , $v_{l,j}:=v(t_{l},x_{j})$ and $w_{l,j}:=w(t_{l},x_{j})$ .

The Crank-Nicolson (CN) finite difference method approximates the spacial derivatives by center difference,

	$\displaystyle\left.v_{x}\right\|_{t=t_{l},x=x_{j}}$	$\displaystyle=$	$\displaystyle\frac{v_{l,j+1}-v_{l,j-1}}{2\Delta x}$
	$\displaystyle\left.v_{xx}\right\|_{t=t_{l},x=x_{j}}$	$\displaystyle=$	$\displaystyle\frac{v_{l,j+1}-2v_{l,j}+v_{l,j-1}}{\Delta x^{2}}$

Along the time dimension, consider a middle point $(t,x_{j}),t\in(t_{l},t_{l+1})$ . The CN method approximates time derivative by center difference and averages the right hand side of equation (13) as

\left.v_{t}\right|_{x=x_{j}}=\frac{v_{l+1,j}-v_{l,j}}{\Delta t}=\phi\left(w_{l% +1,j}+s_{l+1,j}\right)+\psi\left(w_{l,j}+s_{l,j}\right),\psi=1-\phi

CN method effectively sets $\phi=\frac{1}{2}$ . The above difference equation is a general form that has the flavor of choosing $0\leq\phi\leq 1$ from fully explicit to fully implicit. Then differential equation (13) is approximated by the following finite difference equation,

	$\displaystyle\frac{v_{l+1,j}}{\Delta t}-$
	$\displaystyle\phi\left[f_{l+1,j}\frac{v_{l+1,j+1}-2v_{l+1,j}+v_{l+1,j-1}}{% \Delta x^{2}}+g_{l+1,j}\frac{v_{l+1,j+1}-v_{l+1,j-1}}{2\Delta x}+h_{l+1,j}v_{l% +1,j}\right]=$
	$\displaystyle\frac{v_{l,j}}{\Delta t}+\psi\left[f_{l,j}\frac{v_{l,j+1}-2v_{l,j% }+v_{l,j-1}}{\Delta x^{2}}+g_{l,j}\frac{v_{l,j+1}-v_{l,j-1}}{2\Delta x}+h_{l,j% }v_{l,j}\right]+$
	$\displaystyle\phi s_{l+1,j}+\psi s_{l,j}$

At boundary states e.g., $x_{0}$ , this difference equation involves node values $v_{l,-1}$ and $v_{l+1,-1}$ . We can introduce “ghost point” $x_{-1}$ and incorporate with boundary conditions. In the general form of Robin conditions, we can express $v_{l,-1}$ as

v_{l,-1}=\frac{2\alpha_{1}\Delta x}{\beta_{1}}v_{l,0}+v_{l,1}-\frac{2\Delta x}% {\beta_{1}}u_{1}(t_{l})

, thereby eliminate nodal values $v_{l,-1}$ and $v_{l+1,-1}$ . For all states $x_{j}$ to satisfy the difference equations, effectively we solve a linear system of the follow type at each $t_{l}$ time step.

\mathbf{B}_{l+1}v_{l+1}=\mathbf{A}_{l}v_{l}+\mathbf{s}

(14)

The CN method has second order accuracy $\mathcal{O}\left(\Delta t^{2}\right)$ and $\mathcal{O}\left(\Delta x^{2}\right)$ in both time and space dimension. It is unconditionaly stable. Finite difference methods work well in the situation when initial value $v(0,x)$ is not smooth or even non-differentiable. In which case, non-smooth $v(0,x)$ will cause grid values oscillating at the first few time steps but soonly being smoothed out as $t$ departs from 0. To improve accuracy, Mathwrist NPL allows users to specify an initial smoothing time $s>0$ . Internally, we will start with a 4-th order accurate finite difference method using the method of lines to approximate $v(t,x)$ to time $t=s$ and then switch back to the regular CN method.

2.3 Spectral Methods

2.3.1 Introduction

Spectral methods approximate the solution $v(t,x)$ by an expansion function $u(t,x)=\sum_{k=-\infty}^{\infty}a_{k}(t)\phi_{k}(x)$ in terms of trial basis functions $\phi_{k}(x)$ . The basis functions are global orthogonal polynomials. If $u(x)$ is a periodic function, the trial functions are typically choosen to be complex exponentials $\phi_{k}(x)=e^{ikx}$ . If $u(x)$ is not periodic, $\phi_{k}(x)$ are usually choosen to be Chebyshev or Legendre polynomials. With these choices of basis functions, it can be shown that the approximation converges to the true solution at exponential rate, i.e. the coefficient series $a_{k}$ converges to 0 exponentially. This property allows us to use a truncated series approximation.

For simplicity, consider a homogeneous PDE $v_{t}=\mathcal{A}(v)$ and approximate $v(t,x)$ by truncated series expansion,

u^{n}(t,x)=\sum_{k=-n/2}^{n/2}a_{k}(t)\phi_{k}(x)

, where the notation $u^{n}(t,x)$ indicates it is an approximation polynomial of degree $n$ . Different spectral methods apply different techniques to reduce the residual term $\frac{\partial}{\partial t}u^{n}-\mathcal{A}(u^{n})$ . One approach, known as weak formulation, is to select a set of test functions $\psi_{k}(x),k=-\frac{n}{2},\cdots,\frac{n}{2}$ and require

\int_{a}^{b}\left[\frac{\partial}{\partial t}u^{n}+\mathcal{A}(u^{n})\right]% \psi_{k}(x)dx=0,\forall k

An alternative approach takes a set of collocation points $x_{j}\in(a,b)$ and requires the PDE being exactly satisfied at the collocation points,

\left.\frac{\partial}{\partial t}u^{n}\right|_{x=x_{j}}=\mathcal{A}(u^{n}),\forall j

This approach is called spectral collocation method or pseudo-spectral method. It is a strong formulation because it perserves the original differential order. Spectral collocation methods are similar to finite difference methods in the sense that they all impose the strong formulation of the original problem to a set of discrete points.

2.3.2 Chebyshev Collocation

Mathwrist NPL implements a Chebyshev collocation method where we choose Chebyshev polynomials as basis functions

\phi_{k}(x)=T_{k}(x)=\cos\left(k\cos^{-1}x\right)

and collocation points being the extrema points of $n$ -th order Chebyshev polynomial $x_{j}=\cos\frac{\pi j}{n},j=0,\cdots,n$ . The approximation polynomial is

u^{n}(t,x)=\sum_{k=0}^{n}a_{k}(t)T_{k}(x),k=0,\cdots,n

(15)

Coefficients $a_{k}(t)$ can be expressed in terms of node values $u(t,x_{j})$ [5]. We can then impose the strong form of the 1-d parabolic PDE in (10) to $u^{n}(t,x)$ and obtain an ODE form,

\left.\frac{\partial u^{n}}{\partial t}\right|_{x=x_{j}}=f(t,x)\frac{\partial^% {2}u^{n}}{\partial x^{2}}+g(t,x)\frac{\partial u^{n}}{\partial x}+h(t,x)u^{n}+% s(t,x)

(16)

Solving the PDE at collocation points then becomes to solve a system of ODE, where the spatial partial derivatives can be computed from (15) or its equivalent interpolation polynomial. At boundary points $x_{0}$ and $x_{n}$ , nodal values $u_{0}$ and $u_{n}$ need incorporate with boundary conditions. When solving ODE system (16), Mathwrist NPL uses its spectral IDeC method with implicit Euler as the base method.

Not suprisingly, spectral collocation method outperforms when $v(t,x)$ is very smooth in spatial dimension, including initial values $v(0,x)$ . When $v(t,x)$ has fluctuations or strong curvature in states, users can specify the initial smoothing time parameter $s$ like our CN method. In which case, Mathwrist NPL will internally use an higher order finite difference method of lines to approximate $v(t,x)$ over the starting period $0<t<s$ and then switch to spectral collocation method by taking $v(s,x)$ as initial values for the rest of time.

2.4 Method of Lines with Finite Difference

2.4.1 Introduction

The technique of method of lines (MOL) in general sense refers to the family of numerical algorithms that solve a parabolic PDE (10) by solving a system of ODE in the form of (11). Spectral collocation method of course falls into this category. Alternatively, we can approximate the spatial partial derivative operator $\mathcal{A}(v)$ by finite difference. Mathwrist NPL provides 2nd order and 4th order accurate finite difference scheme to approximate $\mathcal{A}(v)$ . Classic IDeC methods are then used to solve the system of ODE with implicit Euler as the base method.

2.4.2 Fourth Order Finite Difference Method

Let $\delta^{(d)}(v)$ be the $d$ -th order finite difference operator wrt spatial state $x$ at interior node $x_{j}$ . Using the method of undetermined coefficients, see [7] pp 19-20, at each grid point $x_{k}$ we approximate $\delta^{(1)}(v)=\frac{\partial v}{\partial x}$ and $\delta^{(2)}(v)=\frac{\partial^{2}v}{\partial x^{2}}$ as linear combinations of grid values $v_{j-2},v_{j-1},v_{j},v_{j+1},v_{j+2}$ .

\displaystyle\left.\delta^{(d)}(v)\right|_{x=x_{j}}

\displaystyle=

\displaystyle c^{(d)}_{0}v_{j-2}+c^{(d)}_{1}v_{j-1}+c^{(d)}_{2}v_{j}+c^{(d)}_{% 3}v_{j+1}+c^{(d)}_{4}v_{j+2}

By matching the coefficients of derivative terms with the 4-th order Taylor expansion around $x_{j}$ , we can obtain

	$\displaystyle c^{(1)}_{0}$	$\displaystyle=$	$\displaystyle\frac{1}{12\Delta x},c^{(1)}_{1}=\frac{-2}{3\Delta x},c^{(1)}_{2}% =0,c^{(1)}_{3}=\frac{2}{3\Delta x},c^{(1)}_{4}=\frac{-1}{12\Delta x}$
	$\displaystyle c^{(2)}_{0}$	$\displaystyle=$	$\displaystyle\frac{-1}{12\Delta x^{2}},c^{(2)}_{1}=\frac{4}{3\Delta x^{2}},c^{% (2)}_{2}=\frac{-5}{2\Delta x^{2}},c^{(2)}_{3}=\frac{4}{3\Delta x^{2}},c^{(2)}_% {4}=\frac{-1}{12\Delta x^{2}}$

This is the 5-point center difference scheme, which approximates the first derivative $v_{x}$ and the second derivative $v_{xx}$ to 4th order accuracy for interior node $x_{j},j=2,\cdots,J-2$ . For near boundary states $x_{1}$ and $x_{J-1}$ , alternative high order finite difference formulas are used. Lastly, high-order accurate schemes are used to approximate $v_{x}$ that appears in Neumann and Robin boundary conditions. Overall the method is 4th order accurate in states and relies on IDeC method to raise the accuracy along the time dimension, i.e. 3 error correction iterations to obtain $\mathcal{O}(h^{4})$ accuracy.

2.4.3 Second Order Finite Difference Method

We could choose the exactly same technique as the CN method in spatial discretization to approximate $\mathcal{A}(v)$ and boundary conditions. Along the time dimension, we then use implicit IDeC to solve the ODE system (11). Compared to the CN method, this approach offers users more freedom to choose spatial discretization $\Delta x$ and time marching step size $h$ . For example, the method affords to use a relatively larger step size $h$ yet retain the same level of accuracy as the CN method.

References

[1] Richard L. Burden, J. Douglas Faires and Annette M. Burden: Numerical Analysis, 10th edition, Cengage Learning, 2016
[2] Pedro E. Zadunaisky: On the Estimation of Errors Propagated in the Numerical Integration of Ordinary Differential Equations, Numer. Math. 27, 21 - 39 (1976), Springer-Verlag 1976
[3] Alok Dutt, Leslie Greengard, Vladimir Rokhlin: Spectral Deferred Correction Methods for Ordinary Differential Equations, BIT, 40(2):241–266, 2000
[4] Anders C. Hansen, John Strain: On the Order of Deferred Correction, Applied Numerical Mathematics, Volume 61, Issue 8, August 2011, Pages 961-973
[5] C. Canuto, M.Y. Hussaini, A. Quarteroni and T. A. Zang: Spectral Methods, Fundamentals in Single Domains, Springer, 2006
[6] John P. Boyd: Chebyshev and Fourier Spectral Methods, Dover, 2nd edition, 2016
[7] J. W. Thomas: Numerical Partial Differential Equations: Finite Difference Methods, Springer, 2010
[8] Sandro Salsa: Partial Differential Equations in Action, from Modelling to Theory, Springer 2008

Differential Equations Mathwrist White Paper Series