MATH383: Numerical methods for differential equations

Overview

• Motivation

• Euler's method

• Runge-Kutta method

• Applications:

  • SIR model

  • Van der Pol oscillator

  • Double pendulum

Motivation

• Many differential equations of practical interest are nonlinear and do not have a closed-form analytical solution.

• Examples:

  Van der Pol oscillator $x^{\text{'}\text{'}}-\mu (1-x^2)x^{\text{'}}+x=0$

  Double pendulum

  $\begin{array}{ll}{\dot{\theta }}_1=\frac{6}{m{l}^2}\frac{2p_1-3\cos ({\theta }_1-{\theta }_2)p_2}{16-9{\cos }^2({\theta }_1-{\theta }_2)}& {\dot{\theta }}_2=\frac{6}{m{l}^2}\frac{8p_1-3\cos ({\theta }_1-{\theta }_2)p_2}{16-9{\cos }^2({\theta }_1-{\theta }_2)}\\ {\dot{p}}_1=-\frac{m{l}^2}{2}({\dot{\theta }}_1{\dot{\theta }}_2\sin ({\theta }_1-{\theta }_2)+3\frac{g}{l}\sin {\theta }_1)& {\dot{p}}_2=\frac{m{l}^2}{2}({\dot{\theta }}_1{\dot{\theta }}_2\sin ({\theta }_1-{\theta }_2)-\frac{g}{l}\sin {\theta }_1)\end{array}$

• Though apparently complicated such systems are of the general form ${𝒚}^{\text{'}}=𝒇(x,𝒚)$

• Methods for systems are essentially the same as those for the scalar DE $y^{\text{'}}=f(x,y)$

• We'll study a few key concepts on numerical solution of $y^{\text{'}}=f(x,y)$

Basic idea

• IVP: $y^{\text{'}}=f(x,y)$, $y\left(x_0\right)=y_0$, find solution over interval $[x_0,x_n]$

• Discretize the interval with step size $h=(x_n-x_0)/n$, $x_k=x_0+kh,k=0,\dots ,n$

• The slope is sampled at points $(x_k,y_k)$, notation $f_k=f(x_k,y_k)$

• The derivative is approximated, e.g., Euler's method $y_k\cong y\left(x_k\right)$, $y_k^{\text{'}}\cong y^{\text{'}}\left(x_k\right)$

  ${\displaystyle y^{\text{'}}\left(x_k\right)=f(x_k,y\left(x_k\right))\cong y_k^{\text{'}}=\frac{y_{k+1}-y_k}{h}=f(x_k,y_k)\Rightarrow }$ ${\displaystyle y_{k+1}=y_k+hf(x_k,y_k)}$

• An approximate solution is built up iteratively

  ${\displaystyle \begin{array}{l}{y}_1=y_0+hf(x_0,y_0)\\ y_2=y_1+hf(x_1,y_1)\\ ⋮\\ y_n=y_n+hf(x_n,y_n)\end{array}}$

Basic analysis of Euler's method: truncation error

• Truncation error. Since $\frac{y_{k+1}-y_k}{h}$ is an approximation of $y_k^{\text{'}}$, an error is introduced at each step of the iteration. The truncation error can be estimated from a Taylor series

  ${\displaystyle \begin{array}{l}{y}_{k+1}=y_k+hf(x_k,y_k)\Rightarrow y\left(x_{k+1}\right)=y\left(x_k\right)+hf(x_k,y_k)+{\tau }_k\\ \mathrm{Error}:e_k=y\left(x_k\right)-y_k=\mathrm{exact}-\mathrm{approximation}(\mathrm{Note}:e_0=0)\\ \mathrm{Since}{x}_{k+1}=x_k+h:\\ y\left(x_{k+1}\right)=y\left(x_k\right)+y^{\text{'}}\left(x_k\right)h+\frac{1}{2}{y}^{\text{'}\text{'}}\left({\xi }_k\right)h^2\\ e_{k+1}=e_k+\frac{1}{2}{y}^{\text{'}\text{'}}\left({\xi }_k\right)h^2\end{array}}$

• Cummulative error (global truncation error)

  ${\displaystyle \begin{array}{l}{e}_1=e_0+\frac{1}{2}{y}^{\text{'}\text{'}}\left({\xi }_0\right)h^2\\ e_2=e_1+\frac{1}{2}{y}^{\text{'}\text{'}}\left({\xi }_1\right)h^2\\ ⋮\\ e_n=e_{n-1}+\frac{1}{2}{y}^{\text{'}\text{'}}\left({\xi }_{n-1}\right)h^2\end{array}\Rightarrow e_n=\sum _{i=0}^{n-1}\frac{1}{2}{y}^{\text{'}\text{'}}\left({\xi }_i\right)h^2⩽max\left|y^{\text{'}\text{'}}\right|\frac{x_n-x_0}{2}h=Kh}$

Basic analysis of Euler's method: stability

• Numerical methods use floating point numbers, a finite-precision approximation of reals

• A numerical method that amplifies inherent rounding error is unstable, e.g.,$\pi -3.141592=e$

• Consider $y^{\text{'}}=f\left(y\right)\cong \lambda y$ approximated by Euler's method

  ${\displaystyle \begin{array}{l}{y}_{k+1}-y_k=\lambda h{y}_k\Rightarrow y_{k+1}=(1+\lambda h)y_k\\ y\left(x_{k+1}\right)=y\left(x_k\right)+\lambda hy\left(x_k\right)+\frac{1}{2}{y}^{\text{'}\text{'}}\left({\xi }_k\right)h^2\end{array}\Rightarrow e_{k+1}=(1+\lambda h)e_k+\frac{1}{2}{y}^{\text{'}\text{'}}\left({\xi }_k\right)h^2}$

  $$

• Euler's method error

  ${\displaystyle \begin{array}{l}{e}_1=(1+\lambda h)e_0+\frac{1}{2}{y}^{\text{'}\text{'}}\left({\xi }_0\right)h^2\\ e_2=(1+\lambda h)e_1+\frac{1}{2}{y}^{\text{'}\text{'}}\left({\xi }_1\right)h^2\\ ⋮\\ e_n=(1+\lambda h)e_{n-1}+\frac{1}{2}{y}^{\text{'}\text{'}}\left({\xi }_{n-1}\right)h^2\end{array}\Rightarrow e_n={(1+\lambda h)}^{n-1}{e}_1+\dots }$

• If $|1+\lambda h|>1$ errors get amplified, Euler's method is unstable, $|1+\lambda h|⩽1$ stable

Runge-Kutta methods

• Euler's method approximates the average slope over interval $[x_k,x_{k+1}]$ by $f(x_k,y_k)$, i.e., at the start of the interval. $y\left(x_{k+1}\right)=y\left(x_k\right)+y^{\text{'}}\left({\xi }_k\right)h$

• Runge-Kutta idea: approximate the average slope by a weighted average

  ${\displaystyle \begin{array}{l}{k}_1=hf(x_i,y_i)\\ k_2=hf(x_i+\frac{h}{2},y_i+\frac{1}{2}{k}_1)\\ k_3=hf(x_i+\frac{h}{2},y_i+\frac{1}{2}{k}_2)\\ k_4=hf(x_i+h,y_i+k_3)\\ y_{i+1}=y_i+\frac{1}{6}(k_1+2k_2+2k_3+k_4)\end{array}}$

• The above method is one particular example of a large class of Runge-Kutta algorithms

• The (global) truncation error is $e_n=K{h}^4$, much more accurate than that for Euler's method ($e_n=Lh$) since $h$ can be considered small (subunitary)

Application 1: SIR model

• Recall the SIR epidemic model: $S^{\text{'}}=-\beta IS$, $I^{\text{'}}=\beta IS-\gamma I$, $R^{\text{'}}=\gamma I$

(%i2) 

plotdf([-beta*I*S,beta*I*S-gamma*I],[S,I],
[trajectory_at,1000,1],
[S,1,1000],[I,1,1000],
[direction,forward],
[parameters,"beta=0.0001,gamma=0.01"],[sliders,"beta=.00001:.01,gamma=0.01:.1"],[versus_t,1])$

(%i3) 

Application 2: Van der Pol oscillator

• $x^{\text{'}\text{'}}-\mu (1-x^2)x^{\text{'}}+x=0\Rightarrow$

  ${\displaystyle x^{\text{'}}=y,y^{\text{'}}=\mu (1-x^2)y-x}$

(%i4) 

plotdf([y,mu*(1-x^2)*y-x],[x,y],
[trajectory_at,1,1],
[x,-5,5],[y,-5,5],
[direction,forward],
[parameters,"mu=1"],[sliders,"mu=-1:1"],[versus_t,1])$

(%i5) 

Application 3: Double pendulum

• $\begin{array}{ll}{\dot{\theta }}_1=\frac{6}{m{l}^2}\frac{2p_1-3\cos ({\theta }_1-{\theta }_2)p_2}{16-9{\cos }^2({\theta }_1-{\theta }_2)}& {\dot{\theta }}_2=\frac{6}{m{l}^2}\frac{8p_1-3\cos ({\theta }_1-{\theta }_2)p_2}{16-9{\cos }^2({\theta }_1-{\theta }_2)}\end{array}$

• Choose $m{l}^2=6$

(%i6) 

plotdf([
 (2*p1-3*cos(th1-th2)*p2)/
 (16-9*(cos(th1-th2))^2),
 (8*p1-3*cos(th1-th2)*p2)/
 (16-9*(cos(th1-th2))^2)],[th1,th2],
[trajectory_at,1,1],
[th1,-4*%pi,4*%pi],
[th2,-4*%pi,4*%pi],
[direction,forward],
[parameters,"p1=1,p2=1"],[sliders,"p1=-1:1,p2=-1:1"],[versus_t,1])$

(%i7)