MATH566 Lesson 23: Numerical ODE

Overview

• ODE review

• First-order ODE initial value problem

• Differential operator approximation from polynomial interpolants

  • forward Euler scheme

  • backward Euler scheme

  • Leapfrog scheme

• Schemes based upon numerical quadrature

  • Adams-Bashforth schemes

  • Adams-Moulton schemes

• Forward Euler analysis

• General analysis techniques: convergence = consistency + stability

ODE review

• An $n^{\mathrm{th}}$-order ordinary differential equation given in explicit form

  ${\displaystyle y^{\left(n\right)}=\frac{d^{n}y}{d{t}^n}=f(t,y,y^{\text{'}},\dots ,y^{(n-1)}).}$

• The problem is to determine the function $y:ℝ\to ℝ$, $y\in C^{\left(n\right)}\left(ℝ\right)$

• $y$ is not given directly but as an equality between two operators acting on $y$

  ${\displaystyle ℒy=f(t,y,y^{\text{'}},\dots ,y^{(n-1)}),ℒ=\frac{d^n}{d{t}^n},f:{ℝ}^{n+1}\to ℝ.}$

• An $n^{\mathrm{th}}$ order ODE is equivalent to a system of $n$ first-order ODEs

  ${\displaystyle {𝒛}^{\text{'}}=𝑭(t,𝒛)}$

  where

  ${\displaystyle 𝒛={\left[\begin{array}{lllll}{z}_1& z_2& \dots & z_{n-1}& z_n\end{array}\right]}^T={\left[\begin{array}{lllll}y& y^{\text{'}}& \dots & y^{(n-2)}& y^{(n-1)}\end{array}\right]}^T,}$ ${\displaystyle 𝑭(t,𝒛)={\left[\begin{array}{lllll}{z}_2\left(t\right)& z_3\left(t\right)& \dots & z_n\left(t\right)& f(t,z_1\left(t\right),\dots ,z_n\left(t\right))\end{array}\right]}^T.}$

  This leads to the central role of numerical solution of first-order ODEs.

First-order ODE initial value problem

• The initial value problem (IVP) for $y:ℝ\to ℝ$ is

  ${\displaystyle y^{\text{'}}=f(t,y),y\left(0\right)=y_0.}$

• IVP has a unique solution for $(t,y)\in [0,T]\times [y_1,y_2]$ if $f$ Lipschitz-continuous

  ${\displaystyle \exists K\in {ℝ}_{+}\mathrm{such}\mathrm{that}|f(t,y_2)-f(t,y_1)|⩽K|y_2-y_1|.}$

  $|f(t,y_2)-f(t,y_1)|=𝒪\left(|y_2-y_1|\right)$ $|f_2-f_1|$ goes to zero as $|y_2-y_1|$.

• Lipschitz is more restrictive than continuity ($\forall \epsilon >0,\exists {\delta }_{\epsilon }$ such that $|t_2-t_1|<{\delta }_{\epsilon }\Rightarrow |y\left(t_2\right)-y\left(t_1\right)|<\epsilon$ (no relation between how |$t_2-t_1$|,$|y_2-y_1|$ go to zero

Figure 1. To solve an ODE IVP is to find a specific integral curve.

Numerical ODE - derivative approximation

• Approaches to numerical solution of $ℒy=f(t,y)$, $ℒ=d/dt$

  • Approximation of the differentiation operator $ℒ$, $ℒ\cong 𝒟$

  • Approximation of the nonlinear operator $f$, $f\cong g$

  • Approximation of the equality between effect of two operators $𝒟y=g(t,y)$

• Methods from approximating $ℒ=d/dt$, $f_i=f(t_i,y_i)$

  • Forward Euler (an explicit method, next $y$ value given directly)

    ${\displaystyle \frac{dy}{dt}\left(t_i\right)\cong \frac{y\left(t_{i+1}\right)-y\left(t_i\right)}{t_{i+1}-t_i}=\frac{y_{i+1}-y_i}{h}\Rightarrow y_{i+1}=y_i+h{f}_i,i=0,1,2,\dots }$

  • Backward Euler (an implicit method, must solve an equation to find $y_i$)

    ${\displaystyle \frac{dy}{dt}\left(t_i\right)\cong \frac{y\left(t_i\right)-y\left(t_{i-1}\right)}{t_i-t_{i-1}}=\frac{y_i-y_{i-1}}{h}\Rightarrow y_i=y_{i-1}+hf(t_i,y_i),i=1,2,\dots }$

  • Leapfrog (centered finite difference)

    ${\displaystyle \frac{dy}{dt}\left(t_i\right)\cong \frac{y_{i+1/2}-y_{i-1/2}}{t_{i+1/2}-t_{i-1/2}}=f_i\Rightarrow y_{i+1/2}=y_{i-1/2}+h{f}_i,}$

Numerical ODE - Adams-Bashforth

• Integrate $y^{\text{'}}=f(t,y)$ over time step $[t_i,t_{i+1}]$, $y_{i+1}-y_i={\int }_{t_i}^{t_{i+1}}f(t,y\left(t\right))dt$

• Approximate $f$ on data set $𝒟=\{(t_i,y_i),(t_{i-1},y_{i-1}),\dots ,(t_{i+1-s},y_{i+1-s})\}$

  ${\displaystyle f(t,y\left(t\right))\cong \sum _{k=1}^s{\ell }_k\left(t\right)f_k,f_k=f(t_{i+1-k},y\left(t_{i+1-k}\right))\cong f(t_{i+1-k},y_{i+1-k}).}$

• The resulting schemes are known as Adams-Bashforth explicit methods

  ${\displaystyle y_{i+1}=y_i+\sum _{k=1}^s(\underset{t_i}{\overset{t_{i+1}}{\int }}{\ell }_k\left(t\right)dt)f_k=y_i+h\sum _{k=1}^s{b}_k{f}_{i+1-k},}$

  $\begin{array}{lllll}s& b_1& b_2& b_3& b_4\\ 1& 1& & & \\ 2& & -& & \\ 3& & -& & \end{array}$

  Table 1. Adams-Bashforth scheme coefficients.

Numerical ODE - Adams-Moulton

• Integrate $y^{\text{'}}=f(t,y)$ over time step $[t_i,t_{i+1}]$, $y_{i+1}-y_i={\int }_{t_i}^{t_{i+1}}f(t,y\left(t\right))dt$

• Approximate $f$ on data set $𝒟=\{(t_{i+1},y_{i+1}),(t_i,y_i),\dots ,(t_{i+2-s},y_{i+2-s})\}$

  ${\displaystyle f(t,y\left(t\right))\cong \sum _{k=0}^{s-1}{\ell }_k\left(t\right)f_k,f_k=f(t_{i+1-k},y\left(t_{i+1-k}\right))\cong f(t_{i+1-k},y_{i+1-k}).}$

• The resulting schemes are known as Adams-Moulton implicit methods

  ${\displaystyle y_{i+1}=y_i+\sum _{k=1}^s(\underset{t_i}{\overset{t_{i+1}}{\int }}{\ell }_k\left(t\right)dt)f_k=y_i+h\sum _{k=1}^s{b}_k{f}_{i+1-k},}$

Numerical analysis of forward Euler (FE)

• Introduce error at step $i$, $e_i=y\left(t_i\right)-y_i.$ ($y\left(t_i\right)$ denotes the exact value)

  ${\displaystyle e_{i+1}-e_i=y\left(t_{i+1}\right)-y\left(t_i\right)-(y_{i+1}-y_i)\Rightarrow }$ ${\displaystyle e_{i+1}-e_i=h{y}^{\text{'}}\left(t_i\right)+\frac{h^2}{2}{y}^{\text{'}\text{'}}\left({\xi }_i\right)-y\left(t_i\right)-hf(t_i,y\left(t_i\right))=\frac{h^2}{2}{y}^{\text{'}\text{'}}\left({\xi }_i\right).}$

• At each step forward Euler introduces an error (the one-step error)

  ${\displaystyle {\tau }_i=e_{i+1}-e_i=\frac{h^2}{2}{y}^{\text{'}\text{'}}\left({\xi }_i\right).}$

• After $N$ steps

  ${\displaystyle e_N-e_0=\frac{h^2}{2}\sum _{i=1}^N{y}^{\text{'}\text{'}}\left({\xi }_i\right).}$

• Exact start $e_0=0$. Obtain for $T=Nh$, that FE is first-order of accuracy.

  ${\displaystyle e_N⩽\frac{N{h}^2}{2}{||y^{\text{'}\text{'}}||}_{\infty }=h\frac{T}{2}{||y^{\text{'}\text{'}}||}_{\infty }=𝒪\left(h\right)}$

Numerical analysis of backward Euler (BE)

• Introduce error at step $i$, $e_i=y\left(t_i\right)-y_i.$ ($y\left(t_i\right)$ denotes the exact value)

  ${\displaystyle e_{i+1}-e_i=y\left(t_{i+1}\right)-y\left(t_i\right)-(y_{i+1}-y_i)\Rightarrow }$ ${\displaystyle e_{i+1}-e_i=y\left(t_{i+1}\right)-y\left(t_{i+1}\right)+h{y}^{\text{'}}\left(t_{i+1}\right)-\frac{h^2}{2}{y}^{\text{'}\text{'}}\left({\xi }_i\right)-hf(t_{i+1},y\left(t_{i+1}\right)).}$

• At each step backward Euler introduces an error (the one-step error)

  ${\displaystyle {\tau }_i=e_{i+1}-e_i=-\frac{h^2}{2}{y}^{\text{'}\text{'}}\left({\xi }_i\right).}$

• After $N$ steps

  ${\displaystyle e_N-e_0=\frac{h^2}{2}\sum _{i=1}^N{y}^{\text{'}\text{'}}\left({\xi }_i\right).}$

• Exact start $e_0=0$. Obtain for $T=Nh$, that BE is first-order of accuracy.

  ${\displaystyle e_N⩽\frac{N{h}^2}{2}{||y^{\text{'}\text{'}}||}_{\infty }=h\frac{T}{2}{||y^{\text{'}\text{'}}||}_{\infty }=𝒪\left(h\right)}$

Numerical analysis of leapfrog (LF)

• Introduce error at step $i$, $e_i=y\left(t_i\right)-y_i.$ ($y\left(t_i\right)$ denotes the exact value)

  ${\displaystyle e_{i+1}-e_i=y\left(t_{i+1}\right)-y\left(t_i\right)-(y_{i+1}-y_i)\Rightarrow }$ ${\displaystyle {\tau }_i=e_{i+1}-e_i=y\left(t_{i+1}\right)-y\left(t_i\right)-hf(t_{i+1/2},y\left(t_{i+1/2}\right))\Rightarrow }$ ${\displaystyle {\tau }_i=y\left(t_{i+1/2}\right)+\frac{h}{2}{y}^{\text{'}}\left(t_{i+1/2}\right)+\frac{h^2}{4}{y}^{\text{'}\text{'}}\left(t_{i+1/2}\right)+\frac{h^3}{8}{y}^{\text{'}\text{'}\text{'}}\left({\xi }_i\right)-\dots }$ ${\displaystyle \dots y\left(t_{i+1/2}\right)+\frac{h}{2}{y}^{\text{'}}\left(t_{i+1/2}\right)-\frac{h^2}{4}{y}^{\text{'}\text{'}}\left(t_{i+1/2}\right)+\frac{h^3}{8}{y}^{\text{'}\text{'}\text{'}}\left({\eta }_i\right)-hf(t_{i+1/2},y\left(t_{i+1/2}\right))}$

• At each step leapfrog introduces an error (the one-step error)

  ${\displaystyle {\tau }_i=e_{i+1}-e_i=\frac{h^3}{4}[\frac{y^{\text{'}\text{'}\text{'}}\left({\xi }_i\right)+y^{\text{'}\text{'}\text{'}}\left({\eta }_i\right)}{2}]=\frac{h^3}{8}{y}^{\text{'}\text{'}\text{'}}\left({\zeta }_i\right).}$

• After $N$ steps with $T=Nh$, that LF is second-order of accuracy.

  ${\displaystyle e_N⩽\frac{N{h}^3}{8}{||y^{\text{'}\text{'}\text{'}}||}_{\infty }=h^2\frac{T}{8}{||y^{\text{'}\text{'}}||}_{\infty }=𝒪\left(h^2\right)}$

Accuracy is not the whole story!

• Consider a small initial error in forward Euler $\stackrel{\sim }{y_0}=y_0+ϵ$ (e.g., floating point)

  ${\displaystyle y^{\text{'}}=\lambda y\Rightarrow y_{i+1}=y_i+h\lambda y_i=(1+h\lambda )y_i=(1+z)y_i}$

• After $N$ steps FE with a perturbed initial condition gives

  ${\displaystyle {\stackrel{\sim }{y}}_N={(1+z)}^N(y_0+ϵ)=y_N+{(1+z)}^Nϵ=y_N+e_N}$

• For $|1+z|>1$ the error $e_N$ increases without bound. The forward Euler scheme is said to be unstable.

• How to approach this? First, precisely define a convergent sequence of approximations ${\left\{y_n\right\}}_{n\in \mathbb{N}}$ for the solution $y\left(t_n\right)$ of the IVP $y^{\text{'}}=f(t,y),y\left(0\right)=y_0$, over the time interval $[0,T]$, $t_n=nh$, $h=T/N$. A numerical method (scheme) is said to be convergent if

  ${\displaystyle {lim}_{\begin{array}{l}h\to 0\\ Nh=T\end{array}}{y}_N=y\left(T\right).}$ ${\displaystyle }$

• Analytical calculations of the above limit are however difficult.

An alternative characterization of convergence

• Consider the model problem $y^{\text{'}}=\lambda y,y\left(0\right)=y_0,\lambda ⩽0$ with solution

  ${\displaystyle y\left(t\right)=e^{\lambda t}{y}_0\Rightarrow y\left(t_n\right)=e^{n\lambda h}{y}_0^{}.}$

• Now, consider a perturbation of the initial conditions

  ${\displaystyle \stackrel{\sim }{y}\left(T\right)=e^{\lambda T}(y_0+\delta )\Rightarrow \epsilon =\stackrel{\sim }{y}\left(T\right)-y\left(T\right)=e^{\lambda T}\delta .}$

• Error $\epsilon$ is maintained small if $\lambda ⩽0$. How does a numerical scheme behave?

• Forward Euler: $y_N={(1+z)}^N{y}_0$. Exponential decay of analytical solution only if

  ${\displaystyle -\frac{2}{\lambda }>h>0.}$

• The above is known as a stability condition.

• In the limit of $h\to 0$, the error $e_N⩽h\frac{T}{2}{||y^{\text{'}\text{'}}||}_{\infty }$ also goes to zero. This is known as consistency.

• In general a numerical scheme is convergent iff it is stable and consistent.