Problem 4 - Numerical Calculus

The precision and stability of numerical methods for differential equation are investigated by using

$$\frac{dy}{dx} = \lambda y \; (\hbox{Re}\; \lambda < 0), \quad y(0) = 1$$

as an example problem. For each of the three methods below, compute errors in the first step and find a condition (stability condition) so that the following holds:

$$\lim_{n \rightarrow \infty} \vert y_n\vert = 0$$

where $y_n$ is a numerical solution for $y(nh)$.

1.

Euler method with step-size $h$.

2.

Backward Euler method with step-size $h$.

3.

Trapezoidal method with step-size $h$.

1.

For the Euler method ( $y_{n+1} = y_n + hf(x_n, y_n)$), the error in the first step is:

$y(x_1) - y_1$	$=$	$y(x_0 + h) - \left( y_0 + h \lambda y_0 \right)$
	$=$	$y(x_0) + h \lambda y(x_0) + \frac{h^2}{2} \lambda^2 y(x_0) - y_0 - h \lambda y_0 + O(h^3)$
	$\approx$	$\frac{h^2}{2} \lambda^2$

The stability condition is determined by:

$\vert y_{n+1} \vert = \vert y_n(1 + h \lambda) \vert = \vert y_0(1 + h \lambda)^{n+1} \vert \xrightarrow[n\rightarrow \infty]{} 0$

Which is true if:

$\vert 1 + h \lambda \vert < 1$	$\Leftrightarrow$	$(1 + h Re(\lambda))^2 + h^2 Im(\lambda)^2 < 1$
	$\Leftrightarrow$	$(\frac{1}{h}+Re(\lambda))^2 + Im(\lambda)^2 < \frac{1}{h^2}$

Which means $\lambda \in \mathcal{D}_{][}((-\frac{1}{h},0),\frac{1}{h})$

2.

For the Euler backward method ( $y_{n+1} = y_n + hf(x_{n+1}, y_{n+1})$), the error in the first step is:

$y(x_1) - y_1$	$=$	$y(x_0 + h) - \left( y_0 + h \lambda y_1 \right)$
	$=$	$y(x_0) + h \lambda y(x_0) + \frac{h^2}{2} \lambda^2 y(x_0) - y_0 - h \lambda y_1 + O(h^3)$
	$\approx$	$h \lambda + \frac{h^2}{2} \lambda^2 -h \lambda (y_0 + h \lambda y_0)$
	$\approx$	$-\frac{h^2}{2} \lambda^2$

The stability condition is determined by:

$\vert y_{n+1} \vert \xrightarrow[n\rightarrow \infty]{} 0$

Since:

$y_{n+1} = y_n + h \lambda y_{n+1} \Rightarrow y_{n+1} = y_n \left(\frac{1}{1 - h \lambda}\right) = y_0 \left( \frac{1}{1 - h \lambda}\right)^{n+1}$

It is satisfied for:

$\left\vert \frac{1}{1 - h \lambda} \right\vert < 1$	$\Leftrightarrow$	$\vert 1 - h \lambda \vert > 1$
	$\Leftrightarrow$	$(\frac{1}{h} - Re(\lambda))^2 + Im(\lambda)^2 > \frac{1}{h^2}$

Which means $\lambda \in \mathbb{C}^2 - \mathcal{D}_{[]}((\frac{1}{h},0),\frac{1}{h})$.

3.

For the trapezoidal method ( $y_{n+1} = y_n + h \frac{f(x_n,y_n) + f(x_{n+1},y_{n+1})}{2}$), the error in the first step is:

$y(x_1) - y_1$	$=$	$y(x_0 + h) - y_0 - h \lambda \frac{y_0 + y_1}{2}$
	$=$	$y(x_0) + h \lambda y(x_0) + \frac{h^2}{2} \lambda^2 y(x_0) + \frac{h^3}{3!}\lambda^3 y(x_0) - y_0 - h \lambda \frac{y_0 + y_1}{2} + 0(h^4)$
	$\approx$	$h \lambda + \frac{h^2}{2} \lambda^2 + \frac{h^3}{3!} \lambda^3 - h \lambda \frac{y_0 + y_0 + h \lambda y_0 + \frac{h^2}{2} \lambda^2 y_0}{2}$
	$\approx$	$\frac{h^3}{3!} \lambda^3 - \frac{h^3}{4} \lambda^3$
	$\approx$	$-\frac{h^3}{12} \lambda^3$

The stability condition is determined by:

$\vert y_{n+1} \vert \xrightarrow[n\rightarrow \infty]{} 0$

Since:

$y_{n+1} = y_n + h \lambda \frac{y_n + y_{n+1}}{2} \Leftrightarrow y_{n+1} = y_... ... y_0 \left( \frac{1 + \frac{h\lambda}{2}}{1 - \frac{h\lambda}{2}} \right)^{n+1}$

It is satisfied for:

$\left\vert \frac{1 + \frac{h\lambda}{2}}{1 - \frac{h\lambda}{2}} \right\vert < 1$	$\Leftrightarrow$	$\left( 1 + \frac{h\lambda}{2} \right)^2 < \left( 1 - \frac{h\lambda}{2} \right)^2$
	$\Leftrightarrow$	$(1+\frac{hRe(\lambda)}{2})^2 + (\frac{hIm(\lambda)}{2})^2 < (1-\frac{hRe(\lambda)}{2})^2 + (\frac{hIm(\lambda)}{2})^2$
	$\Leftrightarrow$	$(1+\frac{hRe(\lambda)}{2})^2 < (1-\frac{hRe(\lambda)}{2})^2$

Which is always true since $Re(\lambda) < 0$ and $h \approx 0^+$.