4

I am looking at the following exercise:

Consider the initial value problem

$\left\{\begin{matrix} x''(t)=x(t)\\ x(0)=a\\ x'(0)=b \end{matrix}\right.$

Write it as a system of First-Order ODES with suitable initial values and show that Euler method can get unstable for a great step $(h)$.

That is the solution that the assistant of the prof gave us: $$y_1=x \Rightarrow y_1'=x'=y_2 | y_1(0)=x(0)=a \\ y_2=x' \Rightarrow y_2'=x''=y_1 | y_2(0)=x'(0)=b $$

$$\binom{y_1}{y_2}'=\begin{pmatrix} 0 & 1\\ 1 & 0 \end{pmatrix} \binom{y_1}{y_2} \text{ with } \binom{y_1(0)}{y_2(0)}=\binom{a}{b}$$

Euler method:

$$\binom{y_1^{n+1}}{y_2^{n+1}}=\binom{y_1^n}{y_2^n}+h \begin{pmatrix} 0 & 1\\ 1 & 0 \end{pmatrix} \binom{y_1^n}{y_2^n}=\binom{y_1^n+hy_2^n}{y_2^n+hy_1^n} \\ \binom{y_1^{n+1}}{y_2^{n+1}}= \begin{pmatrix} 1 & h\\ h & 1 \end{pmatrix} \binom{y_1^n}{y_2^n}$$

$$b=-a$$

The exact solution of the initial value problem is ($t^n=nh$)

$\binom{y_1(t)}{y_2(t)}=e^t \binom{a}{-a} \Rightarrow \binom{y_1(t^n)}{y_2(t^n)}=e^{-t^n} \binom{a}{-a}=e^{-nh} \binom{a}{-a}$, $h$ constant, $n \to +\infty, y \to 0$.

$y=e^{-t}$

Euler method

$\binom{y_1^1}{y_2^1}=\begin{pmatrix} 1 & h\\ h & 1 \end{pmatrix} \binom{y_1^0}{y_2^0}=\begin{pmatrix} 1 & h\\ h & 1 \end{pmatrix} \binom{a}{-a}=(1-h) \binom{a}{-a}$

$\binom{y_1^2}{y_2^2}=\begin{pmatrix} 1 & h\\ h & 1 \end{pmatrix} \binom{y_1^1}{y_2^1}=(1-h)^2 \binom{a}{-a}$

$\dots \dots \dots$

$\binom{y_1^n}{y_2^n}=(1-h)^n \binom{a}{-a}$

$|1-h|>1 \Rightarrow h>2$.

First of all, how do we find that $\binom{y_1^1}{y_2^1}=\begin{pmatrix} 1 & h\\ h & 1 \end{pmatrix} \binom{y_1^0}{y_2^0}=\begin{pmatrix} 1 & h\\ h & 1 \end{pmatrix} \binom{a}{-a}=(1-h) \binom{a}{-a},\binom{y_1^2}{y_2^2}=\begin{pmatrix} 1 & h\\ h & 1 \end{pmatrix} \binom{y_1^1}{y_2^1}=(1-h)^2 \binom{a}{-a}$, $\dots $,$\binom{y_1^n}{y_2^n}=(1-h)^n \binom{a}{-a}$?

I found the following:

$\binom{y_1^1}{y_2^1}=\begin{pmatrix} 1 & h\\ h & 1 \end{pmatrix} \binom{a}{b}, \binom{y_1^2}{y_2^2}=\begin{pmatrix} 1+h^2 & 2h\\ 2h & h^2+1 \end{pmatrix} \binom{a}{b}, \binom{y_1^3}{y_2^3}=\begin{pmatrix} 1+3h^2 & h(h^2+3)\\ h(h^2+3) & 1+3h^2 \end{pmatrix} \binom{a}{b}$

Am I wrong? If not, how could we find the general formula?

Also, don't we find the real solution of the system of the First-Order ODES as follows?

$\begin{vmatrix} -\lambda & 1\\ 1& -\lambda \end{vmatrix}=0 \Rightarrow \lambda^2=1 \Rightarrow \lambda=\pm 1$.

Now we are looking for the eigenvectors.

For $\lambda=1$:

$\begin{pmatrix} -1 & 1 \\ 1 & -1 \end{pmatrix} \binom{u}{w}=\binom{0}{0} \Rightarrow \left\{\begin{matrix} -u+w=0\\ u-w=0 \end{matrix}\right. \Rightarrow \left\{\begin{matrix} w=u\\ u=w \end{matrix}\right. \overset{\text{ we set } u=1}{\Rightarrow } u=w=1$

For $\lambda=-1$:

$\begin{pmatrix} 1 & 1 \\ 1 & 1 \end{pmatrix} \binom{u}{w}=\binom{0}{0} \Rightarrow u+w=0 \Rightarrow u=-w \overset{\text{ we set w=-1}}{\Rightarrow } u=1=-w$

So the solution is of the form:

$\binom{y_1}{y_2}=c_1 \binom{1}{1} e^{t}+ c_2 \binom{-1}{1} e^{-t}$

Using the initial conditions, I got the following:

$\binom{y_1}{y_2}= \frac{a+b}{2} \binom{1}{1} e^t+ \frac{b-a}{2} \binom{-1}{1}e^{-t}$

How could we show that for a great step $h$ the method can get unstable?

EDIT: We have $A=\begin {pmatrix} 1&h\\h&1 \end {pmatrix}$ and we write it as $A=SDS^{-1}$ where $D$ is diagonal. Then $A^n=SD^nS^{-1}$.

$$S=\begin {pmatrix} -1&1\\1&1 \end {pmatrix}$$

$$S^{-1}=\begin {pmatrix} \frac{-1}{2}&\frac{1}{2}\\\frac{1}{2}&\frac{1}{2} \end {pmatrix}$$

$$D=\begin {pmatrix} 1-h&0\\0&1+h \end {pmatrix} \Rightarrow D^n=\begin {pmatrix} (1-h)^n&0\\0&(1+h)^n \end {pmatrix}$$

How do we conclude that $ \binom{y_1^n}{y_2^n}=(1-h)^n \binom{a}{b} ?$

EDIT 2: I found that: $$\binom{y_1^n}{y_2^n}=\begin{pmatrix} (1-h)^n \left( \frac{a-b}{2} \right )+(1+h)^n \left( \frac{a+b}{2} \right )\\ \\ (1-h)^n\left( \frac{b-a}{2} \right )+(1+h)^n \left( \frac{a+b}{2} \right ) \end{pmatrix}$$

I have to show that Euler method gets unstable for a great step $h$.

But why does this hold? $y_1^n$ and $y_2^n$ are unbounded for each $h>0$ since $(1+h)^n \to +\infty$, right? Or am I wrong?

evinda
  • 7,823
  • What are the recent changes? The question doesn't seem to have been edited since it was posted. – robjohn Apr 30 '15 at 19:34
  • Why do you say $b=-a$? Under the exact solution, it should be $t=nh$ and you have lost the $e^{-t}$ solution. – Ross Millikan Apr 30 '15 at 19:39
  • any numerical method will be unstable for this problem. the reason is you have a solution $e^{t}\pmatrix{1\-1}$ no matter what the initial condition once you have a component in the direction $\pmatrix{1\-1},$ is going to grow exponentially. – abel Apr 30 '15 at 19:44
  • @RossMillikan The assistant of the prof set it... He did it elsewhise than I did it... – evinda Apr 30 '15 at 20:00

3 Answers3

3

Given $A=\begin {pmatrix} 1&h\\h&1 \end {pmatrix}$, the way to find $A^n$ is to write $A=SDS^{-1}$ where $D$ is diagonal. Then $A^n=SD^nS^{-1}$. Alpha shows $S=\begin {pmatrix} -1&1\\1&1 \end {pmatrix}, S^{-1}=\begin {pmatrix} -1/2&1/2\\1/2&1/2 \end {pmatrix}, D=\begin {pmatrix} 1-h&0\\0&1+h \end {pmatrix}$ so your computed solutions will grow or shrink with $(1+h)^n, (1-h)^n$

Ross Millikan
  • 374,822
  • Why will the computed solutions grow or shrink with $(1+h)^n, (1-h)^n$ ? – evinda May 02 '15 at 17:22
  • Could you explain me how we deduce that $\binom{y_1^n}{y_1^n}=(1-h)^n \binom{a}{b}$ ? – evinda May 02 '15 at 19:18
  • Can we be sure that such a D exists? If that is so, why? – mathreadler May 02 '15 at 19:30
  • 1
    @mathreadler: You can always find a similarity transformation to Jordan normal form which is almost diagonal and is easy to raise to powers. – Ross Millikan May 02 '15 at 20:02
  • 1
    @evinda: We can't deduce that. You will have $y_1^n=(1-h)^nc+(1+h)^nd$ for some constants $c,d$ and similar for $y_2^n$ You have a two dimensional vector space of solutions. The specific solution will depend on the initial conditions. – Ross Millikan May 02 '15 at 20:04
  • Substituting $S,S^{-1},D$ I got:

    $$\binom{y_1^n}{y_2^n}=\begin{pmatrix} (1-h)^n \left( \frac{a-b}{2} \right )+(1+h)^n \left(\frac{a+b}{2} \right )\ \ (1-h)^n \left( \frac{b-a}{2} \right) +(1+h)^n \left( \frac{a+b}{2} \right ) \end{pmatrix}$$

    I have to show that Euler method gets unstable for a great step $h$.

    $$$$ But why does this hold? $y_1^n$ and $y_2^n$ are unbounded for each $h>0$ since $(1+h)^n \to +\infty$, right? Or am I wrong? $$$$

     – evinda May 02 '15 at 22:04
  • @RossMillikan Also do we find the real solution from the given initial value problem $$\left{\begin{matrix} x''(t)=x(t)\ x(0)=a\ x'(0)=b \end{matrix}\right.$$

    or from the system $\binom{y_1}{y_2}'=\bigl(\begin{smallmatrix} 0 & 1\ 1 & 0 \end{smallmatrix}\bigr) \binom{y_1}{y_2} \ \binom{y_1(0)}{y_2(0)}=\binom{a}{b}$ ?

     – evinda May 03 '15 at 15:35
  • 1
    You find the real solution from the initial problem. You have the real solution in your original post under "Using the initial conditions" You are right that the real solution goes off to infinity in general because of the $e^t$ term. That is probably the source of $b=-a$: somehow you know that the solution should not blow up and you have to make sure that term does not contribute. If you look at the Euler approximation, even if $b=-a$ so the $(1+h)^n$ term is not there, when $h \gt 2$ the other term will blow up. – Ross Millikan May 03 '15 at 15:43
  • 1
    I think my last comment answers your Edit 2: the $(1+h)^n$ always goes to infinity, but if you have initial conditions that eliminate it $(b=-a)$ the real solution will decay as $e^{-t}$. If $h \lt 2$ your computed solution will also decay, but if $h \gt 2, |1-h| \gt 1$ and $(1-h)^n$ blows up. – Ross Millikan May 03 '15 at 15:53
  • @RossMillikan I didn't find that $(1-h)^n \to +\infty$ for $h>2$ because I thought that this would be satisfied when $1-h>1$. So do we require that $|1-h|>1$? If so, could you explain me why? – evinda May 03 '15 at 16:08
  • Since $h$ is always positive, we never have $1-h \gt 1$ The solution in $(1-h)^n$ tends to decay away nicely, as it should because it is supposed to follow the $e^{-t}$ in the real solution. The problem comes when $h \gt 2$ because then it blows up (oscillating in sign). That is what we mean by the method becoming unstable-a well behaved solution is modeled by something that diverges. – Ross Millikan May 03 '15 at 16:35
  • @RossMillikan I still haven't understood why we look at the absolute value $|1-h|$ in order to determine where $(1-h)^n$ tends to. Could you explain it to me? – evinda May 06 '15 at 20:47
  • Because, for example, $-1.1^n$ does not converge. As long as $|1-h| \lt 1, (1-h)^n$ will converge to zero. – Ross Millikan May 06 '15 at 21:07
  • How do we relate the fact that $-1.1^n$ does not converge? – evinda May 06 '15 at 23:01
  • The correct solution is $ae^{-t}$. If you follow this to some reasonably large time $T$, the step in your model will be $h=T/n$ as there are $n$ stepts to get there. The solution your Euler model follows will be $a(1-h)^n$ If we want to follow up to $T=100$, as long as $n \gt 50$ the solution will fall towards zero. When $n \lt 50$ the solution will bounce around wildly. – Ross Millikan May 07 '15 at 17:00
2

Set $J=\begin{bmatrix}0&1\\1&0\end{bmatrix}$ and $I$ the corresponding identity. Then $J^2=I$ and the fundamental solution of the system is $$ e^{tJ}=\cosh(t)I+\sinh(t)J $$ The transformation matrix from initial vector to step $n$ of the Euler method can be split the same way into even and odd parts with respect to h $$ (I+hJ)^n=\frac12((1+h)^n+(1-h)^n)I+\frac12((1+h)^n-(1-h)^n)J $$ Applying the initial vector $v=\binom{a}{-a}$, which is an eigenvector of $I$, $J$ and thus also $I+hJ$ and $\exp(tJ)$ with eigenvalues obtained by replacing $J$ with $-1$, results in $$ Jv=-v,\quad e^{tJ}v=e^{-t}v,\quad (1+hJ)^nv=(1-h)^nv. $$ The other eigenspace consists of vectors $w=\binom{a}{a}$ with $$ Jw=w,\quad e^{tJ}w=e^tw,\quad(1+hJ)^nw=(1+h)^nw. $$

Note that solutions starting in the eigenspace of $-1$ fall to zero, however, for $h>2$, the numerical solution not only oscillates but also grows to infinity. $h=3$ might not seem a reasonable step size, but $x''=400x$ separated as $x'=20y,\,y'=20x$ gives as stability requirement $20h<2$ or $h<0.1$ which seems more realistic.

Lutz Lehmann
  • 126,666
  • So is that what I did so far wrong? Also, at the beginning why does it hold that $e^{tJ}=\cosh(t)I+ \sinh(t)J$ ?Then why do we calculate $(I+hJ)^n$? And in general could you explain me what you do from this part till the end? – evinda Apr 30 '15 at 17:04
  • You did nothing wrong, the given solution specializes very early for the eigenspace for $λ=-1$, i.e., $c_1=0$ and $c_2=-a$. You stay with the general solution. I wanted to highlight that matrix powers reduce to number powers when using an eigenbasis. The linear ODE system $y'=Ay$ has the solution $y(t)=e^{tA}y(0)$, where the matrix exponential uses the same power series as the scalar exponential. – Lutz Lehmann Apr 30 '15 at 17:14
  • So we have the following, right? $\binom{y_1^1}{y_2^1}=\begin{pmatrix} 1 & h\ h & 1 \end{pmatrix} \binom{a}{b}, \binom{y_1^2}{y_2^2}=\begin{pmatrix} 1+h^2 & 2h\ 2h & h^2+1 \end{pmatrix} \binom{a}{b}, \binom{y_1^3}{y_2^3}=\begin{pmatrix} 1+3h^2 & h(h^2+3)\ h(h^2+3) & 1+3h^2 \end{pmatrix} \binom{a}{b}$ $$$$ How can we find the general formula of the solution $\left( \frac{y_1^i}{y_2^i}\right)$ that Euler method gives? – evinda Apr 30 '15 at 17:39
  • You will find that, as I wrote in the second big formula, $$\begin{bmatrix}y^n_1\y^n_2\end{bmatrix}=\begin{bmatrix} e_n(h) & o_n(h) \ o_n(h) & e_n (h) \end{bmatrix}\begin{bmatrix}a\b\end{bmatrix}$$ where $e_n(h),,o_n(h)$ are the even and odd parts of $(1+h)^n$. – Lutz Lehmann Apr 30 '15 at 20:46
1

Euler's Method using the relation $$ \begin{bmatrix} u\\v \end{bmatrix}' = \begin{bmatrix} 0&1\\1&0 \end{bmatrix} \begin{bmatrix} u\\v \end{bmatrix} $$ would give $$ \begin{bmatrix} u\left((k+1)\frac tn\right)\\v\left((k+1)\frac tn\right) \end{bmatrix} = \begin{bmatrix} u\left(k\frac tn\right)\\v\left(k\frac tn\right) \end{bmatrix} + \frac tn \begin{bmatrix} 0&1\\1&0 \end{bmatrix} \begin{bmatrix} u\left(k\frac tn\right)\\v\left(k\frac tn\right) \end{bmatrix} $$ or $$ \begin{align} \begin{bmatrix} u(t)\\v(t) \end{bmatrix} &= \left( \begin{bmatrix} 1&0\\0&1 \end{bmatrix} + \frac tn \begin{bmatrix} 0&1\\1&0 \end{bmatrix} \right)^n \begin{bmatrix} u(0)\\v(0) \end{bmatrix}\\[6pt] &= \begin{bmatrix} -1&1\\1&1 \end{bmatrix} \left( \begin{bmatrix} 1&0\\0&1 \end{bmatrix} + \frac tn \begin{bmatrix} -1&0\\0&1 \end{bmatrix} \right)^n \begin{bmatrix} -1&1\\1&1 \end{bmatrix}^{-1} \begin{bmatrix} u(0)\\v(0) \end{bmatrix}\\[6pt] &= \begin{bmatrix} -1&1\\1&1 \end{bmatrix} \begin{bmatrix} \left(1-\frac tn\right)^n&0\\ 0&\left(1+\frac tn\right)^n \end{bmatrix} \begin{bmatrix} -1&1\\1&1 \end{bmatrix}^{-1} \begin{bmatrix} u(0)\\v(0) \end{bmatrix}\\[6pt] &= \begin{bmatrix} \frac{\left(1+\frac tn\right)^n+\left(1-\frac tn\right)^n}2 &\frac{\left(1+\frac tn\right)^n-\left(1-\frac tn\right)^n}2\\ \frac{\left(1+\frac tn\right)^n-\left(1-\frac tn\right)^n}2 &\frac{\left(1+\frac tn\right)^n+\left(1-\frac tn\right)^n}2 \end{bmatrix} \begin{bmatrix} u(0)\\v(0) \end{bmatrix}\\[6pt] &\to \begin{bmatrix} \cosh(t)&\sinh(t)\\\sinh(t)&\cosh(t) \end{bmatrix} \begin{bmatrix} u(0)\\v(0) \end{bmatrix}\\ \end{align} $$

robjohn
  • 345,667