<!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">"""Solve u'=-a*u, u(0)=I, for t in (0,T] with steps of dt."""</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 — -->
Theme by <a href="https://github.com/orderedlist">orderedlist</a></small></p>
</footer>
</div>
<script src="js/scale.fix.js"></script>
</body>
</html>