<!--
Automatically generated HTML file from DocOnce source
(https://github.com/hplgit/doconce/)
-->
<html>
<head>
<meta http-equiv="Content-Type" content="text/html; charset=utf-8" />
<meta name="generator" content="DocOnce: https://github.com/hplgit/doconce/" />
<title></title>
<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}
div { text-align: justify; text-justify: inter-word; }
</style>
</head>
<!-- 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 -->
<body>
<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 ---------------------- -->
<p>
<center><h4>Mar 15, 2015</h4></center> <!-- date -->
<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 --------------- -->
</body>
</html>