next up previous
Next: 4.1.3 Explicit Methods Up: 4.1 Initial Value Problems Previous: 4.1.1 Euler-Cauchy Algorithm


4.1.2 Stability and Accuracy of Difference Schemes

Let $\mbox{$\bf y$}(t)$ be the exact solution of a DE, and $\mbox{$\bf e$}(t)$ an error: at time $t_{n}$, the algorithm produces $\mbox{$\bf y$}_{n}+\mbox{$\bf e$}_{n}$.
$\Longrightarrow$$t_{n+1}$? For EC, $
\mbox{$\bf y$}_{n+1}+\mbox{$\bf e$}_{n+1}=\mbox{$\bf y$}_{n}+\mbox{$\bf e$}_{n}
+\mbox{$\bf f$}(\mbox{$\bf y$}_{n}+\mbox{$\bf e$}_{n}) \Delta t
$. Generally,
$\displaystyle \mbox{$\bf y$}_{n+1}+\mbox{$\bf e$}_{n+1}$ $\textstyle =$ $\displaystyle T(\mbox{$\bf y$}_{n}+\mbox{$\bf e$}_{n})$  

Expand around the correct solution:
$\displaystyle T(\mbox{$\bf y$}_{n}+\mbox{$\bf e$}_{n})$ $\textstyle \approx$ $\displaystyle T(\mbox{$\bf y$}_{n})+
\left. \frac{dT(\mbox{$\bf y$})}{d\mbox{$\bf y$}} \right\vert _{\mbox{$\bf y$}_{n}} \cdot \mbox{$\bf e$}_{n}$  

or
$\displaystyle \mbox{$\bf e$}_{n+1}$ $\textstyle \approx$ $\displaystyle \left. \frac{dT(\mbox{$\bf y$})}{d\mbox{$\bf y$}} \right\vert _{\...
..._{n}} \cdot \mbox{$\bf e$}_{n}
\equiv \mbox{${\bf G}$} \cdot \mbox{$\bf e$}_{n}$  

The matrix $\mbox{${\bf G}$}$ is called amplification matrix. All its eigenvalues must be within the unit circle:
$\displaystyle \vert g_{i}\vert$ $\textstyle \leq$ $\displaystyle 1 \,,\;\; {\rm for \; all\; \it i}$  





EXAMPLE: EC + Relaxation equation
$\displaystyle T(y_{n})$ $\textstyle \equiv$ $\displaystyle (1-\lambda \Delta t) \, y_{n}$  

$\Longrightarrow$
$\displaystyle \vert 1-\lambda \Delta t\vert$ $\textstyle \leq$ $\displaystyle 1$  

For $\lambda=1$ this condition is met whenever $\Delta t \leq 2$. $\Longrightarrow$Check the previous figure!


EXAMPLE: EC + Harmonic oscillator
$\displaystyle \mbox{$\bf y$}_{n+1}$ $\textstyle =$ $\displaystyle [\mbox{${\bf I}$}+\mbox{${\bf L}$} \Delta t]\cdot \mbox{$\bf y$}_{n} \equiv T(\mbox{$\bf y$}_{n})$  

The amplification matrix is
$\displaystyle \mbox{${\bf G}$}$ $\textstyle \equiv$ $\displaystyle \left. \frac{dT(\mbox{$\bf y$})}{d\mbox{$\bf y$}} \right\vert _{\mbox{$\bf y$}_{n}} =
\mbox{${\bf I}$}+\mbox{${\bf L}$}\Delta t$  

with eigenvalues $g_{1,2}=1\pm i \omega_{0} \Delta t$, so that
$\displaystyle \vert g_{1,2}\vert$ $\textstyle =$ $\displaystyle \sqrt{1+(\omega_{0} \, \Delta t)^{2}} > 1 \;\;\;{\rm always!!}$  

$\Longrightarrow$ EC applied to the harmonic oscillator is never stable.

next up previous
Next: 4.1.3 Explicit Methods Up: 4.1 Initial Value Problems Previous: 4.1.1 Euler-Cauchy Algorithm
Franz J. Vesely Oct 2005
See also:
"Computational Physics - An Introduction," Kluwer-Plenum 2001