Automatically generated HTML file from DocOnce source (https://github.com/hplgit/doconce/)
<meta http-equiv="Content-Type" content="text/html; charset=utf-8" />
<meta name="generator" content="DocOnce: https://github.com/hplgit/doconce/" />
<meta name="keywords" content="keyword,model problem,exponential decay,mesh in time,$\theta$-rule,numerical scheme,finite difference scheme,numerical experiments,Backward Euler method,Crank-Nicolson method,Forward Euler method,error vs time step">


<style type="text/css">
/* blueish style */

/* Color definitions:  http://www.december.com/html/spec/color0.html
   CSS examples:       http://www.w3schools.com/css/css_examples.asp */

body {
  margin-top: 1.0em;
  background-color: #ffffff;
  font-family: Helvetica, Arial, FreeSans, san-serif;
  color: #000000;
h1 { font-size: 1.8em; color: #1e36ce; }
h2 { font-size: 1.6em; color: #1e36ce; }
h3 { font-size: 1.4em; color: #1e36ce; }
a { color: #1e36ce; text-decoration:none; }
tt { font-family: "Courier New", Courier; }
/* pre style removed because it will interfer with pygments */
p { text-indent: 0px; }
hr { border: 0; width: 80%; border-bottom: 1px solid #aaa}
p.caption { width: 80%; font-style: normal; text-align: left; }
hr.figure { border: 0; width: 80%; border-bottom: 1px solid #aaa}
.alert-text-small   { font-size: 80%;  }
.alert-text-large   { font-size: 130%; }
.alert-text-normal  { font-size: 90%;  }
.alert {
  padding:8px 35px 8px 14px; margin-bottom:18px;
  text-shadow:0 1px 0 rgba(255,255,255,0.5);
  border:1px solid #bababa;
  border-radius: 4px;
  -webkit-border-radius: 4px;
  -moz-border-radius: 4px;
  color: #555;
  background-color: #f8f8f8;
  background-position: 10px 5px;
  background-repeat: no-repeat;
  background-size: 38px;
  padding-left: 55px;
  width: 75%;
.alert-block {padding-top:14px; padding-bottom:14px}
.alert-block > p, .alert-block > ul {margin-bottom:1em}
.alert li {margin-top: 1em}
.alert-block p+p {margin-top:5px}
.alert-notice { background-image: url(https://cdn.rawgit.com/hplgit/doconce/master/bundled/html_images/small_gray_notice.png); }
.alert-summary  { background-image:url(https://cdn.rawgit.com/hplgit/doconce/master/bundled/html_images/small_gray_summary.png); }
.alert-warning { background-image: url(https://cdn.rawgit.com/hplgit/doconce/master/bundled/html_images/small_gray_warning.png); }
.alert-question {background-image:url(https://cdn.rawgit.com/hplgit/doconce/master/bundled/html_images/small_gray_question.png); }

div { text-align: justify; text-justify: inter-word; }



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

<h1 id="table_of_contents">
Table of contents</h2>

<a href="#math:problem"> Mathematical problem </a><br>
<a href="#numerical:problem"> Numerical solution method </a><br>
<a href="#___sec2"> Implementation </a><br>
<a href="#___sec3"> Numerical experiments </a><br>
<a href="#___sec4"> Error vs $latex \Delta t$ </a><br>
<a href="#___sec5"> Bibliography </a><br>

<h1 id="math:problem">Mathematical problem</h1>

We address the initial-value problem


u'(t) = -au(t), \quad t \in (0,T],  $

u(0)  = I,                         

where $latex a$, $latex I$, and $latex T$ are prescribed parameters, and $latex 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.

<h1 id="numerical:problem">Numerical solution method</h1>

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

The $latex \theta$-rule <a href="#Iserles_2009">[1]</a> is used to solve <b>(REF to equation ode not supported)</b> numerically:


u^{n+1} = \frac{1 - (1-\theta) a\Delta t}{1 + \theta a\Delta t}u^n,

for $latex n=0,1,\ldots,N_t-1$. This scheme corresponds to

  <li> The <a href="http://en.wikipedia.org/wiki/Forward_Euler_method" target="_self">Forward Euler</a>     scheme when $latex \theta=0$</li>
  <li> The <a href="http://en.wikipedia.org/wiki/Backward_Euler_method" target="_self">Backward Euler</a>     scheme when $latex \theta=1$</li>
  <li> The <a href="http://en.wikipedia.org/wiki/Crank-Nicolson" target="_self">Crank-Nicolson</a>     scheme when $latex \theta=1/2$</li>

<h1 id="___sec2">Implementation </h1>

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


<div class="highlight" style="background: #f8f8f8"><pre style="line-height: 125%"><span></span><span style="color: #008000; font-weight: bold">def</span> <span style="color: #0000FF">solver</span>(I, a, T, dt, theta):
    <span style="color: #BA2121; font-style: italic">&quot;&quot;&quot;Solve u&#39;=-a*u, u(0)=I, for t in (0,T] with steps of dt.&quot;&quot;&quot;</span>
    dt <span style="color: #666666">=</span> <span style="color: #008000">float</span>(dt)            <span style="color: #408080; font-style: italic"># avoid integer division</span>
    Nt <span style="color: #666666">=</span> <span style="color: #008000">int</span>(<span style="color: #008000">round</span>(T<span style="color: #666666">/</span>dt))     <span style="color: #408080; font-style: italic"># no of time intervals</span>
    T <span style="color: #666666">=</span> Nt<span style="color: #666666">*</span>dt                 <span style="color: #408080; font-style: italic"># adjust T to fit time step dt</span>
    u <span style="color: #666666">=</span> zeros(Nt<span style="color: #666666">+1</span>)           <span style="color: #408080; font-style: italic"># array of u[n] values</span>
    t <span style="color: #666666">=</span> linspace(<span style="color: #666666">0</span>, T, Nt<span style="color: #666666">+1</span>)  <span style="color: #408080; font-style: italic"># time mesh</span>

    u[<span style="color: #666666">0</span>] <span style="color: #666666">=</span> I                  <span style="color: #408080; font-style: italic"># assign initial condition</span>
    <span style="color: #008000; font-weight: bold">for</span> n <span style="color: #AA22FF; font-weight: bold">in</span> <span style="color: #008000">range</span>(<span style="color: #666666">0</span>, Nt):    <span style="color: #408080; font-style: italic"># n=0,1,...,Nt-1</span>
        u[n<span style="color: #666666">+1</span>] <span style="color: #666666">=</span> (<span style="color: #666666">1</span> <span style="color: #666666">-</span> (<span style="color: #666666">1-</span>theta)<span style="color: #666666">*</span>a<span style="color: #666666">*</span>dt)<span style="color: #666666">/</span>(<span style="color: #666666">1</span> <span style="color: #666666">+</span> theta<span style="color: #666666">*</span>dt<span style="color: #666666">*</span>a)<span style="color: #666666">*</span>u[n]
    <span style="color: #008000; font-weight: bold">return</span> u, t

<h1 id="___sec3">Numerical experiments </h1>

A set of numerical experiments has been carried out, where $latex I$, $latex a$, and $latex T$ are fixed, while $latex \Delta t$ and $latex \theta$ are varied. In particular, $latex I=1$, $latex a=2$, $latex \Delta t = 1.25, 0.75, 0.5, 0.1$. Figure <a href="#fig:BE">1</a> contains four plots, corresponding to four decreasing $latex \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 <a href="#fig:CN">2</a> and <a href="#fig:FE">3</a>.

<center>  <div id="fig:BE"></div> 
<hr class="figure">
<center><p class="caption">Figure 1:  The Backward Euler method for decreasing time step values.   </p></center>
<p><img src="https://raw.github.com/hplgit/hplgit.github.com/master/teamods/writing_reports/_static/BE.png" align="bottom" width=800></p>

<center>  <div id="fig:CN"></div> 
<hr class="figure">
<center><p class="caption">Figure 2:  The Crank-Nicolson method for decreasing time step values.   </p></center>
<p><img src="https://raw.github.com/hplgit/hplgit.github.com/master/teamods/writing_reports/_static/CN.png" align="bottom" width=800></p>

<center>  <div id="fig:FE"></div> 
<hr class="figure">
<center><p class="caption">Figure 3:  The Forward Euler method for decreasing time step values.   </p></center>
<p><img src="https://raw.github.com/hplgit/hplgit.github.com/master/teamods/writing_reports/_static/FE.png" align="bottom" width=800></p>

<h1 id="___sec4">Error vs $latex \Delta t$ </h1>

How the error


$latex  E^n = \left(\int_0^T (Ie^{-at} - u^n)^2dt\right)^{\frac{1}{2}} $

varies with $latex \Delta t$ for the three numerical methods is shown in Figure <a href="#fig:error">4</a>.

<div class="alert alert-block alert-warning alert-text-normal">
The data points for the three largest $latex \Delta t$ values in the Forward Euler method are not relevant as the solution behaves non-physically.

<center>  <div id="fig:error"></div> 
<hr class="figure">
<center><p class="caption">Figure 4:  Variation of the error with the time step.   </p></center>
<p><img src="https://raw.github.com/hplgit/hplgit.github.com/master/teamods/writing_reports/_static/error.png" align="bottom" width=400></p>

The $latex E$ numbers corresponding to Figure <a href="#fig:error">4</a> are given in the table below.

<table border="1">
<tr><th align="center">$latex \Delta t$</th> <th align="center">$latex \theta=0$</th> <th align="center">$latex \theta=0.5$</th> <th align="center">$latex \theta=1$</th> </tr>
<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>
<div class="alert alert-block alert-summary alert-text-normal">

<li> $latex \theta =1$: $latex E\sim \Delta t$ (first-order convergence).</li>
<li> $latex \theta =0.5$: $latex E\sim \Delta t^2$ (second-order convergence).</li>
<li> $latex \theta =1$ is always stable and gives qualitatively corrects results.</li>
<li> $latex \theta =0.5$ never blows up, but may give oscillating solutions    if $latex \Delta t$ is not sufficiently small.</li>
<li> $latex \theta =0$ suffers from fast-growing solution if $latex \Delta t$ is    not small enough, but even below this limit one can have oscillating    solutions (unless $latex \Delta t$ is sufficiently small).</li>

<h1 id="___sec5">Bibliography </h1>

 <li> <div id="Iserles_2009"></div> <b>A. Iserles</b>.     <em>A First Course in the Numerical Analysis of Differential Equations</em>,     second edition,     Cambridge University Press,     2009.</li>
 <li> <div id="Langtangen_2014"></div> <b>H. P. Langtangen</b>.     <em>A Primer on Scientific Programming With Python</em>,     fifth edition,     Springer,     2016.</li>