Numerical Analysis/ODE in vector form Exercises

All of the standard methods for solving ordinary differential equations are intended for first order equations. When you need to solve a higher order differential equation, you first convert it to a system of first order of equations. Then you rewrite as a vector form and solve this ODE using a standard method. On this page we demonstrate how to convert to a system of equations and then apply standard methods in vector form.

Reduction to a first order system

(Based on Reduction of Order and Converting a general higher order equation.)

I want to show how to convert higher order differential equation to a system of the first order differential equation. Any differential equation of order n of the form

$f(t,u,u^{\prime },u^{″},\cdots ,u^{(n-1)})=u^{(n)}$

can be written as a system of n first-order differential equations by defining a new family of unknown functions

$y_i=u^{(i-1)}\text{for}i=1,2,...n.$

The n-dimensional system of first-order coupled differential equations is then

$\begin{array}{ccc}{y}_1& =& u\\ y_2& =& u^{\prime }\\ y_3& =& u^{″}\\ & ⋮& \\ y_n& =& u^{(n-1)}.\end{array}$

Differentiating both sides yields

$\begin{array}{ccccc}{y_1}^{\prime }& =& u^{\prime }& =& y_2\\ {y_2}^{\prime }& =& u^{″}& =& y_3\\ {y_3}^{\prime }& =& u^{‴}& =& y_4\\ & ⋮& \\ y_{n^{\prime }}& =& u^{(n)}& =& f(t,y_1,\cdots ,y_n).\end{array}$

We can express this more compactly in vector form

${𝐲}^{\prime }=𝐟(t,𝐲)$

where $y_{i+1}=f_i(t,𝐲)$ for $i<n$ and $f_n(t,𝐲)$ = $f(t,y_1,y_2,\cdots ,y_n).$

Exercise

Consider the second order differential equation $u^{″}+u=0$ with initial conditions $u(0)=1$ and $u^{\prime }(0)=0$. We will use two steps with step size $h=\frac{\pi }{8}$ and approximate the values of $u(\frac{\pi }{4})$ and $u^{\prime }(\frac{\pi }{4}).$

Since the exact solution is $u(t)=\cos (t)$ we have $u(\frac{\pi }{4})=0.707106781$ and $u^{\prime }(\frac{\pi }{4})=-0.707106781$.

Exercise 1: Convert this second order differential equation to a system of first order equations.

Solution:

We have second order differential equation $u^{″}(t)=-u(t)$ with initial conditions $u(0)=1$ and $u^{\prime }(0)=0$.

Let $y_1=u$ and $y_2=u^{\prime }$. Differentiating $y_1$ and $y_2$ gives

${y_1}^{\prime }=u^{\prime }=y_2$,

${y_2}^{\prime }=u^{″}=-u=-y_1.$

Thus we have a system of first order equation in vector form

$\left[\begin{array}{c}{y_1}^{\prime }\\ {y_2}^{\prime }\end{array}\right]=\left[\begin{array}{c}{y}_2\\ -y_1\end{array}\right]=f(t,\left[\begin{array}{c}{y}_1\\ y_2\end{array}\right])$

Exercise 2: Apply the Euler method twice.

Solution:

By the Euler's method, $y_{n+1}=y_n+hf(t_n,y_n)$.

I. First apply to get y₁

$y_1=y_0+hf(t_0,y_0)$.

Now we convert $y_1$ as a vector form, where $t_0=0,h=\frac{\pi }{8}$ and $\left[\begin{array}{c}{y}_{0,1}\\ y_{o,2}\end{array}\right]=\left[\begin{array}{c}1\\ 0\end{array}\right]$, then we have

$\begin{array}{cc}\left[\begin{array}{c}{y}_{1,1}\\ y_{1,2}\end{array}\right]& =\left[\begin{array}{c}{y}_{0,1}\\ y_{0,2}\end{array}\right]+hf(t_0,\left[\begin{array}{c}{y}_{0,1}\\ y_{0,2}\end{array}\right])\\ & =\left[\begin{array}{c}{y}_{0,1}\\ y_{0,2}\end{array}\right]+h\left[\begin{array}{c}{y}_{0,2}\\ -y_{0,1}\end{array}\right]\\ & =\left[\begin{array}{c}1\\ 0\end{array}\right]+h\left[\begin{array}{c}0\\ -1\end{array}\right]\\ & =\left[\begin{array}{c}1+0\\ 0-h\end{array}\right]\\ & =\left[\begin{array}{c}1\\ -h\end{array}\right]\\ & =\left[\begin{array}{c}1\\ -\frac{\pi }{8}\end{array}\right]\\ & =\left[\begin{array}{c}1\\ -0.392699082\end{array}\right].\end{array}$

2. Second apply to get y₂

$y_2=y_1+hf(t_1,y_1)$.

Similarly, we convert $y_2$ as a vector form, where $t_1=t_0+h=0+h=h$, where $h=\frac{\pi }{8}$ and $\left[\begin{array}{c}{y}_{1,1}\\ y_{1,2}\end{array}\right]=\left[\begin{array}{c}1\\ -h\end{array}\right]$, then we have

$\begin{array}{cc}\left[\begin{array}{c}{y}_{2,1}\\ y_{2,2}\end{array}\right]& =\left[\begin{array}{c}{y}_{1,1}\\ y_{1,2}\end{array}\right]+hf(t_1,\left[\begin{array}{c}{y}_{1,1}\\ y_{1,2}\end{array}\right])\\ & =\left[\begin{array}{c}{y}_{1,1}\\ y_{1,2}\end{array}\right]+h\left[\begin{array}{c}{y}_{2,1}\\ -y_{1,1}\end{array}\right]\\ & =\left[\begin{array}{c}1\\ -h\end{array}\right]+h\left[\begin{array}{c}-h\\ -1\end{array}\right]\\ & =\left[\begin{array}{c}1-h^2\\ -h-h\end{array}\right]\\ & =\left[\begin{array}{c}1-h^2\\ -2h\end{array}\right]\\ & =\left[\begin{array}{c}1-\frac{{\pi }^2}{8^2}\\ -2\frac{\pi }{8}\end{array}\right]\\ & =\left[\begin{array}{c}0.8457874321\\ -0.785398163\end{array}\right].\end{array}$

Exercise 3: Apply the Backward Euler method twice.

Solution:

By the Backward Euler's method,$y_{n+1}=y_n+hf(t_{n+1},y_{n+1}).$

I. First apply to get y₁

$y_1=y_0+hf(t_1,y_1)$.

Now we convert $y_1$ as a vector form, where $t_1=t_0+h=0+h=h$, where $h=\frac{\pi }{8}$ and $\left[\begin{array}{c}{y}_{0,1}\\ y_{o,2}\end{array}\right]=\left[\begin{array}{c}1\\ 0\end{array}\right]$, then we have

$\begin{array}{cc}\left[\begin{array}{c}{y}_{1,1}\\ y_{1,2}\end{array}\right]& =\left[\begin{array}{c}{y}_{0,1}\\ y_{0,2}\end{array}\right]+hf(t_1,\left[\begin{array}{c}{y}_{1,1}\\ y_{1,2}\end{array}\right])\\ & =\left[\begin{array}{c}{y}_{0,1}\\ y_{0,2}\end{array}\right]+h\left[\begin{array}{c}{y}_{1,2}\\ -y_{1,1}\end{array}\right]\\ & =\left[\begin{array}{c}1\\ 0\end{array}\right]+\left[\begin{array}{c}h{y}_{1,2}\\ -h{y}_{1,1}\end{array}\right]\\ & =\left[\begin{array}{c}1+h{y}_{1,2}\\ -h{y}_{1,1}\end{array}\right].\end{array}$

Now we have to solve a system of linear equations

$\{\begin{array}{cc}{y}_{1,1}=1+h{y}_{1,2}& \\ y_{1,2}=-h{y}_{1,1}\end{array}$

That is,

$\{\begin{array}{cc}{y}_{1,1}-h{y}_{1,2}=1& \\ h{y}_{1,1}+y_{1,2}=0\end{array}$

Set up augmented matrix to solve this system,

$\left(\begin{array}{ccc}1& -h& 1\\ h& 1& 0\end{array}\right)\Rightarrow \left(\begin{array}{ccc}1& -h& 1\\ 0& 1+h^2& -h\end{array}\right)$

Thus, the solution is $y_{1,1}=\frac{1}{1+h^2},y_{1,2}=-\frac{h}{1+h^2}.$

Plugging in $h=\frac{\pi }{8}$, then we have

$\left[\begin{array}{c}{y}_{1,1}\\ y_{1,2}\end{array}\right]=\left[\begin{array}{c}0.8457874321\\ -0.785398163\end{array}\right].$

2. Second apply to get y₂

$y_2=y_1+hf(t_2,y_2)$.

Similarly, we convert $y_2$ as a vector form, with $t_2=t_1+h=h+h=2h$, where$h=\frac{\pi }{8}$ and $\left[\begin{array}{c}{y}_{1,1}\\ y_{1,2}\end{array}\right]=\left[\begin{array}{c}\frac{1}{1+h^2}\\ -\frac{h}{1+h^2}\end{array}\right]$, then we have

$\begin{array}{cc}\left[\begin{array}{c}{y}_{2,1}\\ y_{2,2}\end{array}\right]& =\left[\begin{array}{c}{y}_{1,1}\\ y_{1,2}\end{array}\right]+hf(t_2,\left[\begin{array}{c}{y}_{2,1}\\ y_{2,2}\end{array}\right])\\ & =\left[\begin{array}{c}\frac{1}{1+h^2}\\ -\frac{h}{1+h^2}\end{array}\right]+hf(2h,\left[\begin{array}{c}{y}_{2,1}\\ y_{2,2}\end{array}\right])\\ & =\left[\begin{array}{c}\frac{1}{1+h^2}\\ -\frac{h}{1+h^2}\end{array}\right]+h\left[\begin{array}{c}{y}_{2,2}\\ -y_{2,1}\end{array}\right]\\ & =\left[\begin{array}{c}\frac{1}{1+h^2}+h{y}_{2,2}\\ -\frac{h}{1+h^2}-h{y}_{2,1}\end{array}\right].\end{array}$

Now we have to solve a system of linear equations

$\{\begin{array}{cc}{y}_{2,1}=\frac{1}{1+h^2}+h{y}_{2,2}& \\ y_{2,2}=-\frac{h}{1+h^2}-h{y}_{2,1}\end{array}$

That is,

$\{\begin{array}{cc}{y}_{2,1}-h{y}_{2,2}=\frac{1}{1+h^2}& \\ h{y}_{2,1}+y_{2,2}=-\frac{h}{1+h^2}\end{array}$

Set up augmented matrix to solve this system,

$\left(\begin{array}{ccc}1& -h& \frac{1}{1+h^2}\\ h& 1& -\frac{h}{1+h^2}\end{array}\right)\Rightarrow \left(\begin{array}{ccc}1& -h& \frac{1}{1+h^2}\\ 0& 1+h^2& -\frac{2h}{1+h^2}\end{array}\right)$

Thus, the solution is $y_{2,1}=\frac{1-h^2}{(1+h^2{)}^2},y_{2,2}=-\frac{2h}{(1+h^2{)}^2}.$

Plugging in $h=\frac{\pi }{8}$, then we have

$\left[\begin{array}{c}{y}_{2,1}\\ y_{2,2}\end{array}\right]=\left[\begin{array}{c}0.63487705\\ -0.589546795\end{array}\right].$

Exercise 4: Apply the Midpoint method twice.

Solution:

By the Midpoint method, $y_{n+1}=y_n+hf(t_n+\frac{1}{2}h,y_n+\frac{1}{2}hf(t_n,y_n)).$

I. First apply to get y₁

$y_1=y_0+hf(t_0+\frac{1}{2}h,y_0+\frac{1}{2}hf(t_0,y_0)).$

Now we convert $y_1$ as a vector form, where $t_0=0,h=\frac{\pi }{8}$ and $\left[\begin{array}{c}{y}_{0,1}\\ y_{0,2}\end{array}\right]=\left[\begin{array}{c}1\\ 0\end{array}\right]$, then we have

$\begin{array}{cc}\left[\begin{array}{c}{y}_{1,1}\\ y_{1,2}\end{array}\right]& =\left[\begin{array}{c}{y}_{0,1}\\ y_{0,2}\end{array}\right]+hf(t_0+\frac{1}{2}h,\left[\begin{array}{c}{y}_{0,1}\\ y_{0,2}\end{array}\right]+\frac{1}{2}hf(t_0,\left[\begin{array}{c}{y}_{0,1}\\ y_{0,2}\end{array}\right]))\\ & =\left[\begin{array}{c}{y}_{0,1}\\ y_{0,2}\end{array}\right]+hf(t_0+\frac{1}{2}h,\left[\begin{array}{c}{y}_{0,1}\\ y_{0,2}\end{array}\right]+\frac{1}{2}h\left[\begin{array}{c}{y}_{0,2}\\ -y_{0,1}\end{array}\right])\\ & =\left[\begin{array}{c}1\\ 0\end{array}\right]+hf(0+\frac{1}{2}h,\left[\begin{array}{c}1\\ 0\end{array}\right]+\frac{1}{2}h\left[\begin{array}{c}0\\ -1\end{array}\right])\\ & =\left[\begin{array}{c}1\\ 0\end{array}\right]+h\left[\begin{array}{c}-\frac{1}{2}h\\ -1\end{array}\right]\\ & =\left[\begin{array}{c}1-\frac{1}{2}{h}^2\\ -h\end{array}\right].\end{array}$

Plugging in $h=\frac{\pi }{8}$, then we have

$\left[\begin{array}{c}{y}_{1,1}\\ y_{1,2}\end{array}\right]=\left[\begin{array}{c}0.922893716\\ -0.392699082\end{array}\right].$

2. Second apply to get y₂

$y_2=y_1+hf(t_1+\frac{1}{2}h,y_1+\frac{1}{2}hf(t_1,y_1)).$

Now we convert $y_1$ as a vector form, where $t_1=t_0+h=0+h=h,h=\frac{\pi }{8}$ and $\left[\begin{array}{c}{y}_{1,1}\\ y_{1,2}\end{array}\right]=\left[\begin{array}{c}1-\frac{1}{2}{h}^2\\ -h\end{array}\right]$, then we have

$\begin{array}{cc}\left[\begin{array}{c}{y}_{2,1}\\ y_{2,2}\end{array}\right]& =\left[\begin{array}{c}{y}_{1,1}\\ y_{1,2}\end{array}\right]+hf(t_1+\frac{1}{2}h,\left[\begin{array}{c}{y}_{1,1}\\ y_{1,2}\end{array}\right]+\frac{1}{2}hf(t_1,\left[\begin{array}{c}{y}_{1,1}\\ y_{1,2}\end{array}\right]))\\ & =\left[\begin{array}{c}{y}_{1,1}\\ y_{1,2}\end{array}\right]+hf(h+\frac{1}{2}h,\left[\begin{array}{c}{y}_{1,1}\\ y_{1,2}\end{array}\right]+\frac{1}{2}h\left[\begin{array}{c}{y}_{1,2}\\ -y_{1,1}\end{array}\right])\\ & =\left[\begin{array}{c}1-\frac{1}{2}{h}^2\\ -h\end{array}\right]+hf(\frac{3}{2}h,\left[\begin{array}{c}1-\frac{1}{2}{h}^2\\ -h\end{array}\right]+\frac{1}{2}h\left[\begin{array}{c}-h\\ -1+\frac{1}{2}{h}^2\end{array}\right])\\ & =\left[\begin{array}{c}1-\frac{1}{2}{h}^2\\ -h\end{array}\right]+hf(\frac{3}{2}h,\left[\begin{array}{c}1-\frac{1}{2}{h}^2-\frac{1}{2}{h}^2\\ -h-\frac{1}{2}h+\frac{1}{4}{h}^3\end{array}\right])\\ & =\left[\begin{array}{c}1-\frac{1}{2}{h}^2\\ -h\end{array}\right]+hf(\frac{3}{2}h,\left[\begin{array}{c}1-h^2\\ -\frac{2}{3}h+\frac{1}{4}{h}^3\end{array}\right])\\ & =\left[\begin{array}{c}1-\frac{1}{2}{h}^2\\ -h\end{array}\right]+h\left[\begin{array}{c}-\frac{2}{3}h+\frac{1}{4}{h}^3\\ -1+h^2\end{array}\right]\\ & =\left[\begin{array}{c}1-\frac{1}{2}{h}^2\\ -h\end{array}\right]+\left[\begin{array}{c}-\frac{2}{3}{h}^2+\frac{1}{4}{h}^4\\ -h+h^3\end{array}\right]\\ & =\left[\begin{array}{c}1-\frac{1}{2}{h}^2-\frac{2}{3}{h}^2+\frac{1}{4}{h}^4\\ -h-h+h^3\end{array}\right]\\ & =\left[\begin{array}{c}1-2h^2+\frac{1}{4}{h}^4\\ -2h+h^3\end{array}\right].\end{array}$

Plugging in $h=\frac{\pi }{8}$, then we have

$\left[\begin{array}{c}{y}_{2,1}\\ y_{2,2}\end{array}\right]=\left[\begin{array}{c}0.694547552011673\\ -0.7248390292562377\end{array}\right].$

Exercise 5: Using the values from the Midpoint method at t = h in exercise3, apply the Two-step Adams-Bashforth method once.

Solution:

By the Midpoint method, we found

$\left[\begin{array}{c}{y}_{1,1}\\ y_{1,2}\end{array}\right]=\left[\begin{array}{c}1-\frac{1}{2}{h}^2\\ -h\end{array}\right].$

Now apply the Two-step Adams-Bashforth method when $n=0$, then by formula

$y_2=y_1+\frac{3}{2}hf(t_1,y_1)-\frac{1}{2}hf(t_0,y_0).$

Vector form is following by

$\begin{array}{cc}\left[\begin{array}{c}{y}_{2,1}\\ y_{2,2}\end{array}\right]& =\left[\begin{array}{c}{y}_{1,1}\\ y_{1,2}\end{array}\right]+\frac{3}{2}hf(t_1,\left[\begin{array}{c}{y}_{1,1}\\ y_{1,2}\end{array}\right])-\frac{1}{2}hf(t_0,\left[\begin{array}{c}{y}_{0,1}\\ y_{0,2}\end{array}\right])\\ & =\left[\begin{array}{c}1-\frac{1}{2}{h}^2\\ -h\end{array}\right]+\frac{3}{2}hf(t_1,\left[\begin{array}{c}1-\frac{1}{2}{h}^2\\ -h\end{array}\right])-\frac{1}{2}hf(t_0,\left[\begin{array}{c}1\\ 0\end{array}\right])\\ & =\left[\begin{array}{c}1-\frac{1}{2}{h}^2\\ -h\end{array}\right]+\frac{3}{2}h\left[\begin{array}{c}-h\\ -1+\frac{1}{2}{h}^2\end{array}\right]-\frac{1}{2}h\left[\begin{array}{c}0\\ -1\end{array}\right]\\ & =\left[\begin{array}{c}1-2h^2\\ -2h+\frac{3}{4}{h}^3\end{array}\right].\end{array}$

Plugging in $h=\frac{\pi }{8},$

$\left[\begin{array}{c}{y}_{2,1}\\ y_{2,2}\end{array}\right]=\left[\begin{array}{c}0.691574862\\ -0.739978812\end{array}\right].$

Reference

http://en.wikipedia.org/wiki/Ordinary_differential_equation

http://www.math.ohiou.edu/courses/math3600/lecture29.pdf

http://www.ohio.edu/people/mohlenka/20131/4600-5600/hw7.pdf