# Extending numerical Euler method to higher order differential equations

If it was a second order DE I can follow the logic however I get very confused when trying to numerically solve a third order DE. I have to estimate $y(0.1)$ using step size 0.1

$$xy''' + xy'' + x^2y' + y = x,\quad with \quad y(0) = 1, y'(0) = 0, y''(0) = 1$$

where
$$y_0 = y(x_0)$$ $$y_{n +1} = y_n + hf(x_n, y_n)$$ $$x_{n +1} = x_n + h$$

Now I can set $u_1=y, u_2=y', u_3=y''$ which gives $u_1'=y'=u_2$ and $u_2'=y''=u_3$ and $u_3'=y'''$

So I can sub these in to get

$$xu_3'+xu_3+x^2u_2+y=x\tag{*}$$ but I can also get $$xu_3'+xu_2'+x^2u_1'+y=x\tag{**}$$

Would they both lead to the same answers? The first one does seem easier to work with.

And assuming I go with the first one, I still do not know how to set up my iteration because I have no idea how to relate them.

Approximate the value of $y (0.2)$ , where $y(x)$ is the solution of the initial-value problem $$y''+ xy' + y = 0, y(0) = 1, y'(0) = 2$$ Substitute $u = y'$ to obtain the system $y' = u, u ' = −xu − y$ Euler’s scheme then gives $$y_{n +1}= y_n + hu_n$$ $$u_{n +1}= u_n + h( −x_nu_n − y_n)$$

If we use the step-size $h = 0.1$ and $y_0 = 1$, $u_0= 2$, we obtain $$y_1 = y_0 + (0.1)u_0 = 1 + (0.1)2 = 1.2$$ $$u_1 = u_0 + (0.1)( −x 0u 0 − y_0) = 2 + (0.1)( −1) = 1.9$$ $$y_2 = y_1 + (0.1)u_1 = 1.2 + (0.1)(1.9) = 1.39$$ $$u_2 = u_1 + (0.1)( −x_1 u_1 − y_1) = 1.9 − (0.1)1.39 = 1.761.1.9 + (0.1)$$ Thus we obtain $$y(0.2) ≈ 1.39$$ and $$y'(0.2) ≈ 1.761$$

There is a theorem which states the following.

Let $u:I⊂ℝ→ℝ$ be given by \begin{align*}u^{(n)} &= F(t,u,u',…,u^{(n-1)})\\ u^{(k)}(t_0)&=u_k^0, \qquad k=0,…,n-1. \end{align*} This is called an explicit, scalar initial value problem (IVP) of order $n$. Then $u$ is also given as solution to the system of first order: \begin{align*} u'_{n-1}(t)&=F(t,u_0(t),u_1(t),…,u_{n-1}(t)), & u_{n-1}(t_0)&=u_{n-1}^0 \\ u'_{n-2}(t)&=u_{n-1}(t), & u_{n-2}(t_0)&=u_{n-2}^0 \\ \vdots& \\ u'_{0}(t)&=u_{1}(t),& u_{0}(t_0)&=u_{0}^0 \end{align*} using the notation $u_0(t)=u(t)$ and $u_k=\frac{d^k}{dt^k}u$.

What this theorem states, is that you can rewrite a IVP of higher order to a system of ODEs of order one. Written as a matrix system, if $F$ is linear, you will get

$$u' = Au + b.$$

Note, that your ODE is in fact a linear ODE.

The reason why this theorem is important, is that all theories that are developed for order one IVPs, $u'=f(t,u)$, can be applied to such systems of order one. Especially you can use all numerical solvers like the Explicit Euler or RK-methods.

In your case (using your first result $(*)$) that equation is: $$\begin{bmatrix} u_1' \\ u_2'\\u_3' \end{bmatrix} = \begin{bmatrix}0 & 1 & 0\\ 0& 0 & 1 \\ -1/x & -x &-1 \end{bmatrix} \begin{bmatrix}u_1 \\ u_2 \\ u_3 \end{bmatrix} +\begin{bmatrix}0\\0\\x\end{bmatrix}$$

Note, that you have to also include the equation like $u_1' = u_2$ in that system. Also note, that you can't have any $y$'s in there anymore.

That system could be solved by any ODE solver, e.g. Explicit Euler: $$u^{(t+1)} = u^{(t)} + h(Au^{(t)}+b).$$

Is it clear how to set up your iteration from here?

Now for your other case $(**)$. You will get the following system, modulo my calculation errors: $$\begin{bmatrix}1 & 0 & 0\\ 0& 1 & 0 \\ x^2 & x & x \end{bmatrix}\begin{bmatrix} u_1' \\ u_2'\\u_3' \end{bmatrix} = \begin{bmatrix}0 & 1 & 0\\ 0& 0 & 1 \\ -1 & 0 &0 \end{bmatrix} \begin{bmatrix}u_1 \\ u_2 \\ u_3 \end{bmatrix} +\begin{bmatrix}0\\0\\1\end{bmatrix}$$

So this is something like $$Mu'=\tilde{A}u+\tilde{b}.$$ If you would multiply by $M^{-1}$, you get the same system as above.

For numerical methods, the second system has the disadvantage, that you can't directly apply most ODE-Solver to that system. So I don't see any reason to write it in this form.

But we can just try to create an Explicit Euler Method for these kind of matrix systems:

Replace $$u_i' \approx \frac{u_i(t+h) - u(t)}{h} = \frac{u_i^{(t+1)} - u^{t}}{h}$$ yields \begin{align*} u_1^{(t+1)} &= u_1^t + hF_1 \\ u_2^{(t+1)} &= u_2^t + hF_2 \\ x^2u_1^{(t+1)}+xu_2^{(t+1)}+xu_3^{(t+1)} &= x^2u_1^{(t)}+xu_2^{(t)}+xu_3^{(t)} + hF_3 \end{align*} And that translates to: $$Mu^{(t+1)}=Mu^{(t)}+h(\tilde{A}u^{(t)}+\tilde{b}).$$

Again, if you would mutlipy by $M^{-1}$ you would get the Explicit Euler Method for the system $u'=Au+b$. The reason why this works, is that we assumed $F$ to be a linear function, and the difference quotient is also a linear operator in $f$.

• You said that 'Also note, that you can't have any y's in there anymore.' but in the example I added to the original question that was done with a y in there. Also I have to solve this manually if that makes any difference to the best method to use. – user492455 Oct 17 '17 at 13:23
• Yes, with that comment I wanted to emphasize, that you replace $y$ with $u_1$. In your equations $(*)$ and $(**)$ you still used the variable $y$. That is no problem, if you know what you did. – P. Siehr Oct 17 '17 at 13:43
• I am confused about the t you have introduced. So suppose I wanted to do the first iteration for y=0.1, what would t be? – user492455 Oct 17 '17 at 14:05

Indeed you would use the first form $$xu_3'+xu_3+x^2u_2+u_1=x$$ so you have exactly one derivative in the equation.

As you only perform one step, you get $$u_1(h)= u_1(0)+hu_2(0)$$ or $$y(0.1)=y(0)+0.1\cdot y'(0)$$ which can be evaluated only using the intial values without using the equation in any way.