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$$

$$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$.

  • $\begingroup$ 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. $\endgroup$ – user492455 Oct 17 '17 at 13:23
  • $\begingroup$ 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. $\endgroup$ – P. Siehr Oct 17 '17 at 13:43
  • $\begingroup$ 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? $\endgroup$ – 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.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.