#TITLE (actually governed by the filename): Experiments with Schemes for Exponential Decay

By '''Hans Petter Langtangen''' (hpl at simula.no)
==== Jul 2, 2016 ====

''Summary.'' This report investigates the accuracy of three finite difference
schemes for the ordinary differential equation $u'=-au$ with the
aid of numerical experiments. Numerical artifacts are in particular
demonstrated.

__TOC__

<!-- !split -->
== Mathematical problem ==

:\begin{align} u'(t) &= -au(t), \quad t \in (0,T], \\ u(0) &= I, \end{align}
where $a$, $I$, and $T$ are prescribed parameters, and $u(t)$ is
the unknown function to be estimated. This mathematical model
is relevant for physical phenomena featuring exponential decay
in time, e.g., vertical pressure variation in the atmosphere,
cooling of an object, and radioactive decay.

== Numerical solution method ==

We introduce a mesh in time with points $0 = t_0 < t_1 \cdots < t_{N_t}=T$.
For simplicity, we assume constant spacing $\Delta t$ between the
mesh points: $\Delta t = t_{n}-t_{n-1}$, $n=1,\ldots,N_t$. Let
$u^n$ be the numerical approximation to the exact solution at $t_n$.

The $\theta$-rule [1]
is used to solve Equation (ode) numerically:

:$u^{n+1} = \frac{1 - (1-\theta) a\Delta t}{1 + \theta a\Delta t}u^n,$
for $n=0,1,\ldots,N_t-1$. This scheme corresponds to

<ul>
<li> The [http://en.wikipedia.org/wiki/Forward_Euler_method Forward Euler]
scheme when $\theta=0$
<li> The [http://en.wikipedia.org/wiki/Backward_Euler_method Backward Euler]
scheme when $\theta=1$
<li> The [http://en.wikipedia.org/wiki/Crank-Nicolson Crank-Nicolson]
scheme when $\theta=1/2$
</ul>

== Implementation ==

The numerical method is implemented in a Python function
[2] <code>solver</code> (found in the [http://bit.ly/29ayDx3 <code>model.py</code>] Python module file):

<syntaxhighlight lang="python">
def solver(I, a, T, dt, theta):
"""Solve u'=-a*u, u(0)=I, for t in (0,T] with steps of dt."""
dt = float(dt)            # avoid integer division
Nt = int(round(T/dt))     # no of time intervals
T = Nt*dt                 # adjust T to fit time step dt
u = zeros(Nt+1)           # array of u[n] values
t = linspace(0, T, Nt+1)  # time mesh

for n in range(0, Nt):    # n=0,1,...,Nt-1
u[n+1] = (1 - (1-theta)*a*dt)/(1 + theta*dt*a)*u[n]
return u, t
</syntaxhighlight>

== Numerical experiments ==

A set of numerical experiments has been carried out,
where $I$, $a$, and $T$ are fixed, while $\Delta t$ and
$\theta$ are varied. In particular, $I=1$, $a=2$,
$\Delta t = 1.25, 0.75, 0.5, 0.1$.
Figure fig:BE contains four plots, corresponding to
four decreasing $\Delta t$ values. The red dashed line
represent the numerical solution computed by the Backward
Euler scheme, while the blue line is the exact solution.
The corresponding results for the Crank-Nicolson and
Forward Euler methods appear in Figures fig:CN
and fig:FE.

[[File:BE.png|frame|800px|alt=BE.png|The Backward Euler method for decreasing time step values. (fig:BE)]] <!-- not yet uploaded to common.wikimedia.org -->

[[File:CN.png|frame|800px|alt=CN.png|The Crank-Nicolson method for decreasing time step values. (fig:CN)]] <!-- not yet uploaded to common.wikimedia.org -->

[[File:FE.png|frame|800px|alt=FE.png|The Forward Euler method for decreasing time step values. (fig:FE)]] <!-- user: Tsven4, filename: FE.png, timestamp: 2016-06-12T16:40:58Z -->

== Error vs $\Delta t$ ==

How the error

:$E^n = \left(\int_0^T (Ie^{-at} - u^n)^2dt\right)^{\frac{1}{2}}$
varies with $\Delta t$ for the three numerical methods
is shown in Figure fig:error.

{{mbox
| type = warning
| textstyle = font-size: 90%;
| text = '''Observe:''' The data points for the three largest $\Delta t$ values in the
Forward Euler method are not relevant as the solution behaves
non-physically.
}}

[[File:Error.png|frame|400px|alt=Error.png|Variation of the error with the time step. (fig:error)]] <!-- user: Avatar, filename: error.png, timestamp: 2005-11-29T09:58:44Z -->

The $E$ numbers corresponding to Figure fig:error
are given in the table below.

<table border="1">
<tr><th align="center">$\Delta t$</th> <th align="center">$\theta=0$</th> <th align="center">$\theta=0.5$</th> <th align="center">$\theta=1$</th> </tr>
<tbody>
<tr><td align="right">   1.25                     </td> <td align="right">   7.4630                   </td> <td align="right">   0.2161                     </td> <td align="right">   0.2440                   </td> </tr>
<tr><td align="right">   0.75                     </td> <td align="right">   0.6632                   </td> <td align="right">   0.0744                     </td> <td align="right">   0.1875                   </td> </tr>
<tr><td align="right">   0.50                     </td> <td align="right">   0.2797                   </td> <td align="right">   0.0315                     </td> <td align="right">   0.1397                   </td> </tr>
<tr><td align="right">   0.10                     </td> <td align="right">   0.0377                   </td> <td align="right">   0.0012                     </td> <td align="right">   0.0335                   </td> </tr>
</tbody>
</table>

{{mbox
| type = Summary.
| textstyle = font-size: 90%;
| text = '''Summary.'''
<ol>
<li> $\theta =1$: $E\sim \Delta t$ (first-order convergence).
<li> $\theta =0.5$: $E\sim \Delta t^2$ (second-order convergence).
<li> $\theta =1$ is always stable and gives qualitatively corrects results.
<li> $\theta =0.5$ never blows up, but may give oscillating solutions
if $\Delta t$ is not sufficiently small.
<li> $\theta =0$ suffers from fast-growing solution if $\Delta t$ is
not small enough, but even below this limit one can have oscillating
solutions (unless $\Delta t$ is sufficiently small).
</ol>
}}

== Bibliography ==

<ol>
<li> '''A. Iserles'''.
''A First Course in the Numerical Analysis of Differential Equations'',
second edition,
Cambridge University Press,
2009.
<li> '''H. P. Langtangen'''.
''A Primer on Scientific Programming With Python'',
fifth edition,
Springer,
2016.
</ol>