How is the Approximation Related to the IVP, if at all?

The Big Questions about any scheme used to approximate an ODE and its solutions are:

(a)

Is Method Convergent?

(b)

if so, How fast does error $\searrow 0$ as $h\searrow 0$?

(c)

Is method Consistent?

(c)

Is method Stable?

(d)

Is this the most computationally-efficient method to solve the ODE?

We begin to learn these concepts and also learn how to prove these for a variety of different schemes.

First, a useful lemma:

Lemma. For any real $x$

$1+x\le e^x$

for any $x\ge -1\quad\quad 0\le (1+x)^m\le e^{mx}\quad m>0$

Proof. Taylor          $e^x=1+x+\displaystyle \underbrace{\frac{x^2}{2}e}_{\mbox {remainder}}\mbox{ always } > 0$ between 0 and $x.$

i.e.

$$0\le 1+x\le 1+x+\frac12 x^2 e^{\xi}=e^{x}.$$

The second statement trivially follows.

$\Box$

For simplicity, assume that $f(x,y)$ satisfies stronger Lipschitz condition:

$$\vert f(x,y_1)-f(x, y_2)\vert\le L\vert y_1-y_2\vert\left\{ \... ... -\infty < y_1, y_2 <\infty\\ x_0\le x\le b \end{array}\right.$$

$L>0$ (this condition is stronger than necessary $\cdots$ it just simplifies proofs and avoids technicalities that can be mastered after this case is understood).

Remark. Forward Euler (see 0.1.2) is not the best ODE integrator: but it is the simplest.

Theorem.

(Convergence for Forward Euler) Assume $Y(x)$ is solution of IVP (with $f$ satisfying the Lipschitz condition) and

$$\vert\vert Y''(x)\vert\vert _{\infty} \le M \forall x\in [x_0, b].$$

$$\Rightarrow \{y_h(x_n)\vert x_0\le x_n\le b\}$$ obtained by Euler method satisfies

$$\max_{x_0\le x_n\le b} \vert Y(x_n)-y_h(x_n)\vert\le e^{(b-x_0)L}\vert e_0\vert$$

$$+\frac{e^{(b-x_0)L}-1}{L} \tau(h)$$

where $\tau(h)=\displaystyle \frac{h}{2} M$ and $L$ is the Lipschitz constant associated with the IVP and

$e_0=Y_0-y_h(x_0).$

If, in addition, $\vert Y_0-y_h(x_0)\vert\le c_1h$     as $h\to 0$

for some $c_1\ge 0$ (e.g. $Y_0=y_0$ for all $h$) $\Rightarrow \exists\quad B\ge 0$ constant

$$\max_{x_0\le x_n\le b} \vert Y(x_n)-y_h(x_n)\vert\le Bh$$

Remark. How big is $B$? Could be large!!

Proof. Let $e_n=Y(x_n)-y(x_n), \quad n\ge 0$.

let $\tau_n=\displaystyle \frac{h}{2} Y''(\xi_n)\quad 0\le n\le N(h)\equiv$ number of steps depends on $h$ for fixed $[x_0, b]$ interval. Also,

$\qquad\qquad\qquad\qquad\; x_n\le \xi\le x_{n+1}$

estimate: $\max_{0\le n\le N-1}\vert\tau_n\vert\le \tau(h)=\displaystyle \frac{h}{2}M$

Now,

$$Y_{n+1}=Y_n+h f(x_n, Y_n)+h\tau_n$$

$$y_{n+1}=y_n+h f(x_n, y_n)\qquad 0\le n\le N(h)-1$$

subtracting:

$$\therefore e_{n+1}=e_n+h [f(x_n, Y_n)-f(x_n, y_n)]+h\tau_n$$ (12)

$$\vert e_{n+1}\vert\le \vert e_n\vert+hL\underbrace {\vert{Y_n-y_n}\vert}_{e_n}+h\vert\tau_n\vert$$

$$\vert e_{n+1}\vert\le (1+hL) \vert e_n\vert+ h\tau(h)\qquad 0\le n\le N(h)-1$$

Apply (12) recursively:

$$\vert e_n\vert\le (1+hL)^n \vert e_0\vert+\{1+(1+hL)+\cdots (1+hL)^{n-1}\}h\tau (h)$$

Recall: $1+r+r^2+\cdots r^{n-1}=\displaystyle \frac{r^n -1}{r-1}\qquad r\ne 1$.

$$\therefore \vert e_n\vert\le (1+L)^n \vert e_0\vert+\Big[\frac{(1+hL)^n-1}{L}\Big]\tau (h)$$

using $(1+hL)^n \le e^{nhL}=e^{(x_n-x_0)L}\le e^{(b-x_0)}L$

implies main result:

$\therefore$

$$\max\vert Y(x_h)-y_n(x_n)\vert\le e^{(b-x_0)L}\vert e_0\vert+\displaystyle \frac{e^{(b-x_0)L}-1}{L}\tau(h)$$

The

$$\max \vert Y(x_n)-y_n(x)\vert\le Bh$$ (13)

follows trivially from $\nearrow$ above, with

$$B=c_1e^{(b-x_0)L}+\displaystyle \bigg[\frac{e^{(b-x_0)L}-1}{L}\bigg]\frac{M}{2}$$

since by assumption $$\vert Y_0-y_h(x_0)\vert\le c_1h$$     as      $h\rightarrow 0$

$\Box$

Equation (13) gives rate of convergence of method. (Parenthetically, a method that does not converge is useless). But $B$ estimate may be too large. We can sharpen the estimate if the following holds:

Corollary.

Same hypothesis as previous theorem. In addition

$$\frac{\partial f}{\partial y}(x,y)\le 0 \left\{\begin{array}{c} x_0\le x\le b\\ \infty <y<\infty \end{array}\right.$$

$\Rightarrow$ For $h$ sufficiently small

$$\vert Y(x_n)-y(x_n)\vert\le \vert e_0\vert+\displaystyle \frac{h}{2}(x_n-x_0)\max \vert Y''(x)\vert$$

for $x_0\le x_n\le b$.

Proof. Apply Mean Value Theorem to

$$e_{n+1}=e_n+h [f(x_n, Y_n)-f(x_n, y_n)]+h\tau_n$$ (14)

$$e_{n+1}=e_n+h \frac{\partial f}{\partial Y} (x_n, \zeta_n) e_n+ \frac{h^2}{2}Y''(\xi_n)$$ (15)

$$\zeta_n$$between$$y_n(x_n)$$ and$$Y(x_n)$$

Since $y_n$ converges to $Y(x_n)$ on $[x_0, b]\Rightarrow \displaystyle \frac{\partial f(x_n, \zeta_n)}{\partial Y}\rightarrow \frac{\partial f(x, Y(x))}{\partial y}$

and thus bounded in magnitude over $[x_0, b]$. Pick $h_0 > 0$ so that

$$1+h\frac{\partial f(x_n, \zeta_n)}{\partial Y}\ge -1\qquad x_0\le x_n \le b\qquad \forall h\le h_0$$

by assumption $$\frac{\partial f}{\partial Y}<0 \forall h$$. Apply these 2 facts to (15) to get

$$\vert e_{n+1}\vert\le \vert e_n\vert+\displaystyle \frac{h^2}{2} \vert Y''(\xi_n)\vert$$

and then by induction show $\vert e_n\vert\le \vert e_0\vert+\displaystyle \frac{h^2}{2}[ \vert Y''(\xi_0)\vert+\cdots +\vert Y''(\xi_{n-1})\vert]$

which leads to result:

$\Box$

Next we consider the ``stability of solutions''. An IVP may posess stable or unstable solutions, or both. Whether it does depends on $f(x,Y)$, and on the initial data. If the IVP problem is approximated using a numerical scheme, we would like the numerical scheme to have the approximate solution behave as the IVP solutions being studied. However, it is possible that the discretization, may behave in ways that are different from the IVP being approximated...the discretization depends not only on $f$ and on the initial data, but also on how the equation is discretized, on how big the step size is, on how the dependent and independent variables are represented mathematically and on the machine. The study of convergence enabled us to determine whether the approximate solution approached the real solution of the IVP as $h$ got smaller for general choices of initial data and at what rate (see0.1.5). Stability will detemine whether systematic errors (which are usually very small) such as round-off and/or uncertainty in the initial data will make arbitrarily close solution paths separate from each other at a rate greater than linear. For the IVP, this is considered its innate behavior and it is important to know whether what you're approximating has this behavior. For the numerical approximation, we want to know if what we are seeing is due to the innate behavior of the IVP or due to using an inappropriate scheme for the numerical approximation of the IVP.

Remark Numerical instability usually leads to spectacularly bad results, i.e. code crashes. But if we had to rank what's worse, lack of convergence or instability, lack of convergence is actually worse: the reason is that lacking convergence means that we are not solving the IVP we think we're solving but some other IVP! Computationally, also, when instabilities manifest themselves they usually force you to do something about it. But sometimes non-convergent numerical schemes are happy to provide you with all sorts of answers and you'd not suspect anything wrong...well, till you kill someone by solving the wrong equation in the first place.

Stability of IVP

$$\mbox {Take } \left\{\begin{array}{ll} y'=f(x,y)+\delta (x) y(x_0)=Y_0+\varepsilon \end{array}\right.$$ (16)

$$\left\{\begin{array}{ll} Y'=f(x,Y) Y(x_0)=Y_0 \end{array}\right.$$ (17)

$f(x,y)$ continuous and satisfies a Lipschitz Condition (L.C.). Also, assume $\delta (x)$ continuous for all $(x,y)$ interior points to convex set. We show that solution to (16) is controlled by $\varepsilon$ and is unique. Then, we'll use this to infer stability and well-poseness. Note: $\varepsilon$ is a perturbation parameter.

Theorem. (Perturbed Problem): Assume hypothesis on $f(x,y)$ as above and hypothesis on $\delta (x)$ as above. (16) has solution $Y(x;\delta, \varepsilon )$ on $[x_0-\alpha, x_0+\alpha], \alpha >0$, uniformly $\forall \varepsilon$ perturbations and $\delta (x)$ that satisfy

$$\vert\varepsilon \vert\le \varepsilon _0 \qquad \vert\vert\delta \vert\vert _{\infty}\le\varepsilon _0$$

for $\varepsilon _0$ sufficiently small. In addition, if $Y(x)$ is solution of (17), then

$$\displaystyle \max_{\vert x-x_0\vert\le \alpha} \vert Y(x)-Y(x;\d... ...vert\le k[\vert\varepsilon \vert +\alpha \vert\vert\delta \vert\vert _{\infty}]$$

$$k=\displaystyle \frac{1}{1-\alpha L}$$

Proof: exercise. This theorem assumes that $\alpha L$ are sufficiently small. You'll also need

$$1+r+r^2+ \ldots = \frac{1}{1-r}$$

for $r < 1$ $\Box$

Studying the stability of an equation enables us to tell whether

1. [1)]

   Forward Euler (or any other numerical scheme) produces approximate solutions that are close enough to $Y(x)$, the exact solution. Very oftenly we need to determine how close as well.

2. [2)]

   To identify whether equation is ``STIFF'' and/or ``ILL-CONDITIONED'' (this is a topic considered a little later, but for now just think of ``Stiff'' as ``very difficult'' to solve numerically.

For simplicity, consider $\varepsilon$ perturbations (the $\delta (x)\pm 0$ case is a little more involved but enters as per previous theorem).

Take $Y_0\rightarrow Y_0+\varepsilon :$

$$Y'(x;\varepsilon ) = f(x, Y(x;\varepsilon )) \quad\quad x_0-\alpha\le x\le x_0+\alpha$$

$$Y(x_0; \varepsilon ) = Y_0+\varepsilon .$$ (18)

Subtract (18) from (17): let $Z(x)=Y(x, \varepsilon )-Y(x) \Rightarrow Z(x_0, \varepsilon )=\varepsilon$

then

$$Z'(x;\varepsilon )=f(x,Y (x;\varepsilon )) -f(x, Y(x))\approx \displaystyle \frac{\partial f}{\partial Y}(x,Y(x))Z(x;\varepsilon )$$ (19)

Is $Y(x;\varepsilon )$ close to $Y(x)$ as $x\rightarrow \infty$? Maybe, if $\varepsilon$ is small enough and $[x_0-\alpha, x_0 +\alpha]$ small.

We can solve (19)

$$Z(x;\varepsilon ) \approx \varepsilon e^{\int^{x}_{x_0}\displaystyle \frac{df}{dY}(t,Y(t))dt}$$

If $$\frac{\partial f}{\partial Y}(t, Y(t))\le 0$$     $\vert x_0-t\vert\le\alpha \Rightarrow z(x,\varepsilon )$ remains bounded by $\varepsilon$ as $x$ increases.

$$\Rightarrow$$   WE SAY THAT $$(\ref{xx2})$$ is WELL-CONDITIONED !

Example:

$\left\{\begin{array}{lll} Y'=\lambda Y +g(x) & \lambda>0 &\lambda\mbox { constant.}\\ Y(0)=Y_0 \end{array}\right.$

$$\frac{\partial f}{\partial Y}=\lambda$$     and $z(x,\varepsilon )=\varepsilon e^{\lambda x}$ exactly

                 $\therefore$ perturbations grow large as $x$ increases.

$\Rightarrow$ ``ILL-CONDITIONED''         and ``STIFF''        if $\lambda$ LARGE

Example:

$$\left\{\begin{array}{ll} Y'=100Y-101 e^{-x} Y(0)=1 \end{array}\right.$$ (20)

   Solution $$Y(x)=e^{-x}$$

take perturbations

$$Y'=100Y-101e^{-x}\\ Y(0)=1+\varepsilon$$

   Solution: $$Y(x;\varepsilon ) =e^{-x}+\varepsilon e^{100 x}$$

which rapidly departs from true solution $Y(x)$. We say that (20) is ill-conditioned.

For well-conditioned we require that $$\int^{x}_{x_0}\displaystyle \frac{\partial f}{\partial Y} (t,y(t))dt$$ be bounded from above by 0 or a small positive number as $x$ increases $\Rightarrow Z(x;\varepsilon )$ will be bounded by constant times $\varepsilon$.

If $$\frac{\partial f}{\partial Y}\le 0$$ but large $\Rightarrow$ call ODE STIFF

and these cases present problems, numerically.

$\Box$

Stiffness is a qualitative assessment of an ODE or a system of them. In a single ODE stiffness can be assessed as we did above by having some good bounding criteria for $f$, and it is the bounding value that determines how ``stiff'' the ODE is. In a system of ODE's stiffness not only brings into play the size of each $f$ but also the relative size of each of these. That is, in addition to the value of each $f_i$, what also comes into play is the wide discrepancy in the rate of change of the $f_i$'s. If you think of $x$ as a time parameter and can bound the rate of change of each $f_i$ by a constant, say, then if these constants are very disparate we say the ODE system is stiff and it manifests its complication in the existence of a wide span of time scales in the behavior of the solution.

Stability Analysis of Forward Euler Scheme

first we motivate problem with important example:
Example) $\left\{\begin{array}{ll} y'=\alpha y, & \alpha \qquad \mbox{constant}.\\ y(0)=y_0 >0 \end{array}\right.$

It has the exact solution: $y(x)=e^{\alpha x} y_0$. Assume $x>0$.

Using Forward Euler: $y_{n+1}=(1+\alpha h)y_n\rightarrow y_{n+1}=(1+\alpha h)^n y_0$ take $h$ constant for simplicity.

Case (a):	$-1<\alpha h<0$	$\Rightarrow$ solution positive and approaching 0.
Case (b):	$\alpha h<-1$	$\Rightarrow$ sign of solution alternatives as $n$ increases.
Case (c):	$\alpha h>0\Rightarrow$	increases at each step

$\Box$

Definition: In general, for $n=0, 1, 2, \cdots$ a scheme of the form $y_{n+1}=g(y_n, y_{n-1}, \cdots y_0)$ is said to be ``explicit.'' If $y_{n+1}=g(y_{n+1}, y_n, y_n \cdots y_0)$ then scheme is said to be ``implicit.'' Example) Forward Euler is said to be an ``explicit'' scheme because each $y_{n+1}$ can be solved in terms of $y_n$.

Stability Analysis for Forward Euler:

Consider

$$\left\{\begin{array}{ll} z_{n+1}=z_n+h [f(x_n, z_n)+\delta (x_n)] & 0\le n\le N(h)-1 z(0)=y_0 +\varepsilon \end{array}\right.$$ (21)

and

$$\left\{\begin{array}{ll} y_{n+1} =y_n+hf(x_n, y_n) y_0=y(x_0) \end{array}\right.$$ (22)

Look at $\{z_n\}$ and $\{y_n\}$    as $h\rightarrow 0$ and as $n$ increases:

let

$$e_n=z_n-y_n\quad , n\ge 0\Rightarrow \quad e_0=\varepsilon$$

Subtract (22) from (21):

$$e_{n+1}=e_n+h[f(x_n, z_n)-f(x_n, y_n)]+h\delta (x_n)$$

has same form as equation (15) $\therefore$

$$\displaystyle \max_{0\le n\le N(h)} \vert z_n-y_n\vert\le e^{(b-x... ...\vert +\Big[\frac{e^{(b-x_{0})L}-1}{L}\Big] \vert\vert\delta\vert\vert _\infty$$

$\therefore \exists k, k_2$ independent of $h$ with

$$\max_{0\le n\le N(h)} \vert z_n-y_n\vert\le k_1 \vert\varepsilon \vert +k_2 \vert\vert\delta\vert\vert _\infty.$$

Effect of Rounding Errors:

take $\rho_n$ to be the local rounding error and let

$$\tilde{y}_{n+1}=\tilde{y}_n+h f(x_n, \tilde{y}_n)+\rho_n \quad n=0,1\ldots N(h)-1$$ (23)

$$\mbox {let }\quad \left\{\begin{array}{l} \tilde{y}_n \mbox{ are... ...\displaystyle \max_{0\le n\le N(h)-1}\quad \vert\rho_n\vert \end{array}\right.$$

$$Y_{n+1}=Y_n +h f(x_n, y_n)+\displaystyle \frac{h^2}{2} y'' (\xi_n)$$ (24)

Subtract (24) from (23): $\tilde {e}_{n+1}=\tilde{e}_n+h[f(x_n, Y_n)-f(x_n, \tilde{y}_n)]+h\tau_n-\rho_n$

where $\tilde{e}_n=Y(x_n)-\tilde{y}(x_n)\quad and\quad \tau_n\equiv \displaystyle \frac{h}{2}Y''(\xi_n)$

Use same arguments as before, but let $\tau_n -\rho_n/h$ replace $\tau_n$ of before

$$\vert\tilde{e}_n\vert\le e^{(b-x_0)L}\bigg\vert Y_0-\tilde{y}_0\b... ...gg\vert\frac{e^{(b-x_0)L}-1}{L} \bigg\vert\Big[\tau(h)+\frac{\rho(h)}{h}\Big].$$

On a finite precision machine $\rho(h)$ will not decrease as $h\to 0\ldots$ it'll remain finite and approximately constant. Take

$$u\equiv \frac{\rho(h)}{\vert\vert Y\vert\vert _{\infty}}$$ then

$$\vert\tilde {e}_n\vert\le c\Big[\frac{h}{2}\big\vert\big\vert Y''... ...ty}+\frac{u}{h}\big\vert \big\vert Y\big\vert\big\vert _\infty\Big]\equiv E(h)$$

We can find the $h$ which minimizes the error. Call it $h^{\star}$. To find, set $$\frac{dE(h^{\star})}{dh}=0$$ and $h^*$ corresponding to minima.

See Figure 5.

Figure 5: Effect of Rounding Errors on Forward Euler

Asymptotic Error Analysis

Recall that if $B(x,h)$ is a function defined for $x_0\le x \le b$, for sufficiently small $h\Rightarrow$

$B(x,h)={\cal O}(h^p)$	$p>0$ near $\exists$ a constant $c$ such that
$\vert B(x,h)\vert\le ch^p$	$x_0\le x \le b$

Theorem. (Euler Error): Assume $Y(x)=$ solution of ODE and 3 times continuously differentiable. Assume $f_y$ and $f_{yy}$ are continuous and bounded for $x_0\le x\le b, -\infty <y<\infty$. Let the initial value $y_n(x_0)$ satisfy

$$Y_0-y_h(x_0)=\delta_{0}h+{\cal O}(h^2)$$

usually this error is machine precision or zero.

Then the error in Forward Euler's $y_{n+1}=y_n+hf(x_n, y_n)$ satisfies

$$Y(x_n)-y_h(x_n)=D(x_n)h+{\cal O}(h^2)$$

where $\left\{\begin{array}{ll} D'(x)=f_y(x, Y(x)) D(x)+\displaystyle \frac12 Y''(x)\\ D(x_0)=\delta_0 \end{array}\right.$

Proof. before the proof, let's do an example:

Example: $y'=-y$     $y(0)=1$

has solution $Y(x)=e^{-x}$. The $D(x)$ equation is

$\left\{\begin{array}{ll} D'=-D+\displaystyle \frac12e^{-x}\\ D(0)=0 \end{array}\right.$

$\therefore D(x)=\displaystyle \frac12 xe^{-x}$

So the error for Forward Euler is $Y(x_n)-y_h(x_n)\approx\displaystyle \frac{h}{2}x_ne^{-x_{n}}$

$\therefore$ the relative error $$\frac{Y(x_n)-y_h(x_n)}{Y(x_n)}\approx\frac{h}{2}x_n$$

Calculation shows that the error is linearly proportional to $h$.

Remark. We say Euler is an $\lq\lq {\cal O}(h)$ method''

Proof. Use Taylor's

$$Y(x_{n+1})=Y(x_n)+hY'(x_n)+\displaystyle \frac{h^2}{2} Y''(x_n)+\displaystyle \frac{h^3}{6}Y''' (\xi_n)$$ (25)

$$x_n\le\xi_n\le x_{n+1}$$

We have

$$Y'(x)=f(x, Y(x))$$ (26)

and

$$y_{n+1}=y_n+hf(x_n, y_n).$$ (27)

Subtract (27) from (25) and use (26)

$$e_{n+1}=e_n+h\big[f(x_n,Y_n)-f(x_n,y_n)\big]+\frac{h^2}{2}Y''(x_n)+ \frac{h^3}{6}Y(''')(\xi_n).$$ (28)

Continuity in $f(x,y)$ allows us to expand around $Y_n$:

$$f(x_n,y_n)=f(x_n,Y_n)+(y_n-Y_n)f_y(x_n,Y_n)+\displaystyle \frac12(y_n-Y_n)^2 f_{yy}(x_n,\xi)$$

$$\mbox {for }\xi_n \mbox{ between } y_n \mbox{and } Y_{n}.$$

Plug into (28)

$$e_{n+1}= [1+hf_y(x_n, Y_n)]e_n+\frac{h^2}{2}Y''(x_n)+B_n$$

$$B_n=\frac{h^3}{6} Y^{(''')}(\xi_n)-\frac12 hf_{yy}(x_n, \xi_n) e^2_n=0(h^3, he^2_n)$$ (29)

Neglecting $B_n$, then the error

$\left\{\begin{array}{ll} e_{n+1}\approx\underbrace{[1+hf_y(x_n, Y_n)]e_n+\frac... ...elta_0h (\mbox {since }Y_0-y_h(x_0)=\delta_0h+{\cal O}(h^2)) \end{array}\right.$ So $e_n={\cal O}(h)$

and $e_{n+1}=D(x_n)h+{\cal O}(h^2)$

Now, need to show that $g_n$ is principal part of error $e_n$. Let

$$k_n=e_n-g_n.\quad k_0=e_0-g_0= {\cal O}(h^2)$$    by$$\left\{\begin{array}{l} Y_0-y(x_0)=\delta_0h+{\cal O}(h^2)\\ g_0=\delta _oh \end{array}\right.$$

$$k_{n+1} =[1+h fy(x_n, Y_n)] k_n+B_n$$

$$k_{n+1}\le (1+hL)\vert k_n\vert + O(h^3)$$

$$\vert k_n\vert ={\cal O}(h)^2 \therefore$$

$$e_n=g_n+k_n=[h{\cal O}(x_n)+ {\cal O}(h^2)]+O(h^2)$$

$\Box$

Forward Euler is a simple but not always appropriate scheme for solving ODE'S. Let's consider some alternatives: