<!doctype html>
<html>
  <head>
    <meta charset="utf-8">
    <meta http-equiv="X-UA-Compatible" content="chrome=1">
    <title>Solve the world's simplest differential equation</title>

    <link rel="stylesheet" href="css/styles.css">
    <meta name="viewport" content="width=device-width, initial-scale=1, user-scalable=no">
    <!--[if lt IE 9]>
    <script src="//html5shiv.googlecode.com/svn/trunk/html5.js"></script>
    <![endif]-->
  </head>
  <body>
    <div class="wrapper">
      <header>
        <h1>Main Permanent Header</h1>
        <p>Permanent SubHeader</p>

        <!-- picture below the heading on the left -->
	<img src="https://raw.github.com/hplgit/doconce/master/doc/design/fig/discrete_function.png" width=220>

        <!-- *Three* potential link buttons below the picture
        <ul>
         <li><a href="...">Button <strong>1</strong></a></li>
         <li><a href="...">Button <strong>2</strong></a></li>
         <li><a href="...">Button <strong>3</strong></a></li>
        </ul>
        -->
      </header>

      <!-- Here goes the main page --->
      <section>


<!-- tocinfo
{'highest level': 1,
 'sections': [(" Solve the world's simplest differential equation ",
               1,
               None,
               '___sec0'),
              (' Mathematical problem ', 2, None, '___sec1'),
              (' Numerical solution method ', 2, None, '___sec2'),
              (' Implementation ', 2, None, '___sec3'),
              (' Numerical experiments ', 2, None, '___sec4')]}
end of tocinfo -->





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



<!-- ------------------- main content ---------------------- -->

<h1 id="___sec0">Solve the world's simplest differential equation </h1>

<h2 id="___sec1">Mathematical problem </h2>

<p>
This exercise addresses the differential equation 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 constant 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 id="___sec2">Numerical solution method </h2>

<p>
Derive the \( \theta \)-rule scheme for solving \eqref{ode} numerically
with time step \( \Delta t \):

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

Here, \( n=0,1,\ldots,N-1 \).

<p>
<b>Hint.</b>\n
Set up the Forward Euler, Backward Euler, and the Crank-Nicolson (or
Midpoint) schemes first. Then generalize to the \( \theta \)-rule above.



<h2 id="___sec3">Implementation </h2>

<p>
The numerical method is implemented in a Python function
<code>solver</code> (found in the <a href="https://github.com/hplgit/INF5620/tree/gh-pages/src/decay/experiments/decay_mod.py" target="_self"><tt>decay_mod</tt></a> module):

<p>

<!-- code=python (!bc pycod) typeset with pygments style "default" -->
<div class="highlight" style="background: #f8f8f8"><pre style="line-height: 125%"><span style="color: #008000; font-weight: bold">from</span> <span style="color: #0000FF; font-weight: bold">numpy</span> <span style="color: #008000; font-weight: bold">import</span> linspace, zeros

<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>
    N <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> N<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(N<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, N<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>, N):    <span style="color: #408080; font-style: italic"># n=0,1,...,N-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
</pre></div>

<h2 id="___sec4">Numerical experiments </h2>

<p>
Fix the values of where \( I \), \( a \), and \( T \).
Then vary \( \Delta t \) for \( \theta=0,1/2,1 \).
Illustrate that if \( \Delta t \) is not sufficiently small,
\( \theta=0 \) and \( \theta=1/2 \) can give non-physical solutions
(more precisely, oscillating solutions).

<p>
Perform experiments and determine empirically the convergence
rate for \( \theta=0,1/2,1 \).

<!-- ------------------- end of main content --------------- -->


      </section>

      <footer>
        <!-- <p><small>Hosted on GitHub Pages &mdash; -->
        Theme by <a href="https://github.com/orderedlist">orderedlist</a></small></p>
      </footer>

    </div>
    <script src="js/scale.fix.js"></script>

  </body>
</html>