<!-- Use MathJax to render mathematics -->
<script type="text/x-mathjax-config">
MathJax.Hub.Config({
  TeX: {
     equationNumbers: {  autoNumber: "AMS"  },
     extensions: ["AMSmath.js", "AMSsymbols.js", "autobold.js"]
  }
});
</script>
<script type="text/javascript"
 src="http://cdn.mathjax.org/mathjax/latest/MathJax.js?config=TeX-AMS-MML_HTMLorMML">
</script>
</head>

<body bgcolor="white">

<title>Numerical investigations</title>

<center><h1>Experiments with Schemes for Exponential Decay</h1></center>

<center><b>Hans Petter Langtangen</b></center>
<center><b>Simula Research Laboratory</b></center>
<center><b>University of Oslo</b></center>
<center><h4>August 20, 2012</h4></center>

<b>Summary.</b> 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.

<h2>Mathematical problem</h2>

We address the initial-value problem

$$
\begin{align}
u'(t) &= -au(t), \quad t \in (0,T], \label{ode}\\
u(0)  &= I,                         \label{initial:value}
\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.


<h2>Numerical solution method</h2>

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

The \( \theta \)-rule is used to solve \eqref{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-1 \). This scheme corresponds to

<ul>
  <li> The Forward Euler scheme when \( \theta=0 \)
  <li> The Backward Euler scheme when \( \theta=1 \)
  <li> The Crank-Nicolson scheme when \( \theta=1/2 \)
</ul>

<h2>Implementation</h2>

The numerical method is implemented in a Python function:
<pre>
def theta_rule(I, a, T, dt, theta):
    """Solve u'=-a*u, u(0)=I, for t in (0,T] with steps of dt."""
    N = int(round(T/float(dt)))  # no of intervals
    u = zeros(N+1)
    t = linspace(0, T, N+1)

    u[0] = I
    for n in range(0, N):
        u[n+1] = (1 - (1-theta)*a*dt)/(1 + theta*dt*a)*u[n]
    return u, t
</pre>

<h2>Numerical experiments</h2>

We define a set of numerical experiments 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 \).


<h3>The Backward Euler method</h3>
<img src="BE.png" width="800">


<h3>The Crank-Nicolson method</h3>
<img src="CN.png" width="800">


<h3>The Forward Euler method</h3>
<img src="FE.png" width="800">



<h3>Error versus time discretization</h3>

How \( E \) varies with \( \Delta t \) for
\( \theta = 0, 0.5, 1 \) is shown below.
<p>
<img="error.png", width="400">

</body>
</html>