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\left(t,u,u',u'',\ \cdots ,\ u^{(n-1)}\right)=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)}\quad {\text{for}}\quad i=1,2,...n\,.$

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

${\begin{array}{rclcl}y_{1}&=&u\\y_{2}&=&u'\\y_{3}&=&u''\\&\vdots &\\y_{n}&=&u^{(n-1)}.\\\end{array}}$

Differentiating both sides yields

${\begin{array}{rclclcl}y_{1}'&=&u'&=&y_{2}\\y_{2}'&=&u''&=&y_{3}\\y_{3}'&=&u'''&=&y_{4}\\&\vdots &\\y_{n}'&=&u^{(n)}&=&f(t,y_{1},\cdots ,y_{n}).\\\end{array}}$

We can express this more compactly in vector form

$\mathbf {y} '=\mathbf {f} (t,\mathbf {y} )$

where $\ y_{i+1}=f_{i}\left(t,\mathbf {y} \right)$ for $i<n$ and $\ f_{n}\left(t,\mathbf {y} \right)$ = $\ f\left(t,y_{1},y_{2},\cdots ,y_{n}\right)\,.$

Exercise

Consider the second order differential equation $\ u''+u=0$ with initial conditions $\ u{(0)}=1$ and $\ u'{(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'{({\frac {\pi }{4}})}.$

Since the exact solution is $u(t)=\cos(t)$ we have $\ u{({\frac {\pi }{4}})}=0.707106781$ and $u'{({\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'{(0)}=0$.

Let $y_{1}=u$ and $y_{2}=u'$. Differentiating $y_{1}$ and $y_{2}$ gives

$\ y_{1}'=u'=y_{2}$,

$\ y_{2}'=u''=-u=-y_{1}.$

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

$\left[{\begin{array}{c}y_{1}'\\y_{2}'\end{array}}\right]=\left[{\begin{array}{c}y_{2}\\-y_{1}\end{array}}\right]=f\left(t,\left[{\begin{array}{c}y_{1}\\y_{2}\end{array}}\right]\right)$

Exercise 2: Apply the Euler method twice.

Solution:

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

I. First apply to get y₁

$y_{1}=y_{0}+hf\left(t_{0},y_{0}\right)$.

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{aligned}\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\left(t_{0},\left[{\begin{array}{c}y_{0,1}\\y_{0,2}\end{array}}\right]\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{aligned}}$

2. Second apply to get y₂

$y_{2}=y_{1}+hf\left(t_{1},y_{1}\right)$.

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{aligned}\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\left(t_{1},\left[{\begin{array}{c}y_{1,1}\\y_{1,2}\end{array}}\right]\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{aligned}}$

Exercise 3: Apply the Backward Euler method twice.

Solution:

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

I. First apply to get y₁

$\ y_{1}=y_{0}+hf\left(t_{1},y_{1}\right)$.

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{aligned}\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\left(t_{1},\left[{\begin{array}{c}y_{1,1}\\y_{1,2}\end{array}}\right]\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}{hy_{1,2}}\\{-hy_{1,1}}\end{array}}\right]\\&=\left[{\begin{array}{c}{1+hy_{1,2}}\\{-hy_{1,1}}\end{array}}\right]\,.\end{aligned}}$

Now we have to solve a system of linear equations

$\ \left\{{\begin{array}{l l}y_{1,1}=1+hy_{1,2}&\quad \\y_{1,2}=-hy_{1,1}\quad \\\end{array}}\right.$

That is,

$\ \left\{{\begin{array}{l l}y_{1,1}-hy_{1,2}=1&\quad \\hy_{1,1}+y_{1,2}=0\quad \\\end{array}}\right.$

Set up augmented matrix to solve this system,

$\left({\begin{array}{cc|c}1&-h&1\\h&1&0\\\end{array}}\right)\Rightarrow \left({\begin{array}{cc|c}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\left(t_{2},y_{2}\right)$.

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{aligned}\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\left(t_{2},\left[{\begin{array}{c}y_{2,1}\\y_{2,2}\end{array}}\right]\right)\\&=\left[{\begin{array}{c}{\frac {1}{1+h^{2}}}\\{-{\frac {h}{1+h^{2}}}}\end{array}}\right]+hf\left(2h,\left[{\begin{array}{c}y_{2,1}\\y_{2,2}\end{array}}\right]\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}}}+hy_{2,2}}\\{-{\frac {h}{1+h^{2}}}-hy_{2,1}}\end{array}}\right]\,.\end{aligned}}$

Now we have to solve a system of linear equations

$\ \left\{{\begin{array}{l l}y_{2,1}={\frac {1}{1+h^{2}}}+hy_{2,2}&\quad \\y_{2,2}=-{\frac {h}{1+h^{2}}}-hy_{2,1}\quad \\\end{array}}\right.$

That is,

$\ \left\{{\begin{array}{l l}y_{2,1}-hy_{2,2}={\frac {1}{1+h^{2}}}&\quad \\hy_{2,1}+y_{2,2}=-{\frac {h}{1+h^{2}}}\quad \\\end{array}}\right.$

Set up augmented matrix to solve this system,

$\left({\begin{array}{cc|c}1&-h&{\frac {1}{1+h^{2}}}\\h&1&-{\frac {h}{1+h^{2}}}\\\end{array}}\right)\Rightarrow \left({\begin{array}{cc|c}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\left(t_{n}+{\frac {1}{2}}h,y_{n}+{\frac {1}{2}}hf\left(t_{n},y_{n}\right)\right).$

I. First apply to get y₁

$y_{1}=y_{0}+hf\left(t_{0}+{\frac {1}{2}}h,y_{0}+{\frac {1}{2}}hf\left(t_{0},y_{0}\right)\right).$

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{aligned}\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\left(t_{0}+{\frac {1}{2}}h,\left[{\begin{array}{c}y_{0,1}\\y_{0,2}\end{array}}\right]+{\frac {1}{2}}hf\left(t_{0},\left[{\begin{array}{c}y_{0,1}\\y_{0,2}\end{array}}\right]\right)\right)\\&=\left[{\begin{array}{c}y_{0,1}\\y_{0,2}\end{array}}\right]+hf\left(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]\right)\\&=\left[{\begin{array}{c}1\\0\end{array}}\right]+hf\left(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]\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{aligned}}$

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\left(t_{1}+{\frac {1}{2}}h,y_{1}+{\frac {1}{2}}hf\left(t_{1},y_{1}\right)\right).$

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{aligned}\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\left(t_{1}+{\frac {1}{2}}h,\left[{\begin{array}{c}y_{1,1}\\y_{1,2}\end{array}}\right]+{\frac {1}{2}}hf\left(t_{1},\left[{\begin{array}{c}y_{1,1}\\y_{1,2}\end{array}}\right]\right)\right)\\&=\left[{\begin{array}{c}y_{1,1}\\y_{1,2}\end{array}}\right]+hf\left(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]\right)\\&=\left[{\begin{array}{c}{1-{\frac {1}{2}}h^{2}}\\{-h}\end{array}}\right]+hf\left({\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]\right)\\&=\left[{\begin{array}{c}{1-{\frac {1}{2}}h^{2}}\\{-h}\end{array}}\right]+hf\left({\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]\right)\\&=\left[{\begin{array}{c}{1-{\frac {1}{2}}h^{2}}\\{-h}\end{array}}\right]+hf\left({\frac {3}{2}}h,\left[{\begin{array}{c}{1-h^{2}}\\{-{\frac {2}{3}}h+{\frac {1}{4}}h^{3}}\end{array}}\right]\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{aligned}}$

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\left(t_{1},y_{1}\right)-{\frac {1}{2}}hf\left(t_{0},y_{0}\right).$

Vector form is following by

${\begin{aligned}\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\left(t_{1},\left[{\begin{array}{c}y_{1,1}\\y_{1,2}\end{array}}\right]\right)-{\frac {1}{2}}hf\left(t_{0},\left[{\begin{array}{c}y_{0,1}\\y_{0,2}\end{array}}\right]\right)\\&=\left[{\begin{array}{c}{1-{\frac {1}{2}}h^{2}}\\{-h}\end{array}}\right]+{\frac {3}{2}}hf\left(t_{1},\left[{\begin{array}{c}{1-{\frac {1}{2}}h^{2}}\\{-h}\end{array}}\right]\right)-{\frac {1}{2}}hf\left(t_{0},\left[{\begin{array}{c}1\\0\end{array}}\right]\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{aligned}}$

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