boost/odeint
2
0
Fork 0
mirror of https://github.com/boostorg/odeint.git synced 2026-01-27 07:02:11 +00:00
Code Issues Projects Releases Wiki Activity
Files
5b48c56451fd39ccfcc4a1fc9c3aff0fd15c2e77
odeint/doc/boost_numeric_odeint/tutorial/stiff_systems.html
Karsten Ahnert 93c1d4cd12 removing observer from docs
2012-06-26 07:09:51 +02:00

166 lines
18 KiB
HTML
Raw Blame History

<html>
<head>
<meta http-equiv="Content-Type" content="text/html; charset=US-ASCII">
<title>Stiff systems</title>
<link rel="stylesheet" href="../../boostbook.css" type="text/css">
<meta name="generator" content="DocBook XSL Stylesheets V1.75.2">
<link rel="home" href="../../index.html" title="Chapter&#160;1.&#160;Boost.Numeric.Odeint">
<link rel="up" href="../tutorial.html" title="Tutorial">
<link rel="prev" href="chaotic_systems_and_lyapunov_exponents.html" title="Chaotic systems and Lyapunov exponents">
<link rel="next" href="special_topics.html" title="Special topics">
</head>
<body bgcolor="white" text="black" link="#0000FF" vlink="#840084" alink="#0000FF">
<table cellpadding="2" width="100%"><tr><td valign="top"></td></tr></table>
<hr>
<div class="spirit-nav">
<a accesskey="p" href="chaotic_systems_and_lyapunov_exponents.html"><img src="../../images/prev.png" alt="Prev"></a><a accesskey="u" href="../tutorial.html"><img src="../../images/up.png" alt="Up"></a><a accesskey="h" href="../../index.html"><img src="../../images/home.png" alt="Home"></a><a accesskey="n" href="special_topics.html"><img src="../../images/next.png" alt="Next"></a>
</div>
<div class="section">
<div class="titlepage"><div><div><h3 class="title">
<a name="boost_numeric_odeint.tutorial.stiff_systems"></a><a class="link" href="stiff_systems.html" title="Stiff systems">Stiff systems</a>
</h3></div></div></div>
<p>
An important class of ordinary differential equations are so called stiff
system which are characterized by two or more time scales of different order.
Examples of such systems are found in chemical systems where reaction rates
of individual sub-reaction might differ over large ranges, for example:
</p>
<p>
<span class="emphasis"><em>d S<sub>&#8203;1</sub> / dt = - 101 S<sub>&#8203;2</sub> - 100 S<sub>&#8203;1</sub></em></span>
</p>
<p>
<span class="emphasis"><em>d S<sub>&#8203;2</sub> / dt = S<sub>&#8203;1</sub></em></span>
</p>
<p>
To solve stiff systems efficiently using numerics the Jacobian
</p>
<p>
<span class="emphasis"><em>J = d f<sub>&#8203;i</sub> / d x<sub>&#8203;j</sub></em></span>
</p>
<p>
is needed. Here is the definition of the above example
</p>
<p>
</p>
<pre class="programlisting"><span class="keyword">typedef</span> <span class="identifier">boost</span><span class="special">::</span><span class="identifier">numeric</span><span class="special">::</span><span class="identifier">ublas</span><span class="special">::</span><span class="identifier">vector</span><span class="special">&lt;</span> <span class="keyword">double</span> <span class="special">&gt;</span> <span class="identifier">vector_type</span><span class="special">;</span>
<span class="keyword">typedef</span> <span class="identifier">boost</span><span class="special">::</span><span class="identifier">numeric</span><span class="special">::</span><span class="identifier">ublas</span><span class="special">::</span><span class="identifier">matrix</span><span class="special">&lt;</span> <span class="keyword">double</span> <span class="special">&gt;</span> <span class="identifier">matrix_type</span><span class="special">;</span>
<span class="keyword">struct</span> <span class="identifier">stiff_system</span>
<span class="special">{</span>
<span class="keyword">void</span> <span class="keyword">operator</span><span class="special">()(</span> <span class="keyword">const</span> <span class="identifier">vector_type</span> <span class="special">&amp;</span><span class="identifier">x</span> <span class="special">,</span> <span class="identifier">vector_type</span> <span class="special">&amp;</span><span class="identifier">dxdt</span> <span class="special">,</span> <span class="keyword">double</span> <span class="comment">/* t */</span> <span class="special">)</span>
<span class="special">{</span>
<span class="identifier">dxdt</span><span class="special">[</span> <span class="number">0</span> <span class="special">]</span> <span class="special">=</span> <span class="special">-</span><span class="number">101.0</span> <span class="special">*</span> <span class="identifier">x</span><span class="special">[</span> <span class="number">0</span> <span class="special">]</span> <span class="special">-</span> <span class="number">100.0</span> <span class="special">*</span> <span class="identifier">x</span><span class="special">[</span> <span class="number">1</span> <span class="special">];</span>
<span class="identifier">dxdt</span><span class="special">[</span> <span class="number">1</span> <span class="special">]</span> <span class="special">=</span> <span class="identifier">x</span><span class="special">[</span> <span class="number">0</span> <span class="special">];</span>
<span class="special">}</span>
<span class="special">};</span>
<span class="keyword">struct</span> <span class="identifier">stiff_system_jacobi</span>
<span class="special">{</span>
<span class="keyword">void</span> <span class="keyword">operator</span><span class="special">()(</span> <span class="keyword">const</span> <span class="identifier">vector_type</span> <span class="special">&amp;</span> <span class="comment">/* x */</span> <span class="special">,</span> <span class="identifier">matrix_type</span> <span class="special">&amp;</span><span class="identifier">J</span> <span class="special">,</span> <span class="keyword">const</span> <span class="keyword">double</span> <span class="special">&amp;</span> <span class="comment">/* t */</span> <span class="special">,</span> <span class="identifier">vector_type</span> <span class="special">&amp;</span><span class="identifier">dfdt</span> <span class="special">)</span>
<span class="special">{</span>
<span class="identifier">J</span><span class="special">(</span> <span class="number">0</span> <span class="special">,</span> <span class="number">0</span> <span class="special">)</span> <span class="special">=</span> <span class="special">-</span><span class="number">101.0</span><span class="special">;</span>
<span class="identifier">J</span><span class="special">(</span> <span class="number">0</span> <span class="special">,</span> <span class="number">1</span> <span class="special">)</span> <span class="special">=</span> <span class="special">-</span><span class="number">100.0</span><span class="special">;</span>
<span class="identifier">J</span><span class="special">(</span> <span class="number">1</span> <span class="special">,</span> <span class="number">0</span> <span class="special">)</span> <span class="special">=</span> <span class="number">1.0</span><span class="special">;</span>
<span class="identifier">J</span><span class="special">(</span> <span class="number">1</span> <span class="special">,</span> <span class="number">1</span> <span class="special">)</span> <span class="special">=</span> <span class="number">0.0</span><span class="special">;</span>
<span class="identifier">dfdt</span><span class="special">[</span><span class="number">0</span><span class="special">]</span> <span class="special">=</span> <span class="number">0.0</span><span class="special">;</span>
<span class="identifier">dfdt</span><span class="special">[</span><span class="number">1</span><span class="special">]</span> <span class="special">=</span> <span class="number">0.0</span><span class="special">;</span>
<span class="special">}</span>
<span class="special">};</span>
</pre>
<p>
</p>
<p>
The state type has to be a <code class="computeroutput"><span class="identifier">ublas</span><span class="special">::</span><span class="identifier">vector</span></code>
and the matrix type must by a <code class="computeroutput"><span class="identifier">ublas</span><span class="special">::</span><span class="identifier">matrix</span></code>
since the stiff integrator only accepts these types. However, you might want
use non-stiff intgrators on this system, too - we will do so later for demonstration.
Therefore we want to use the same function also with other state_types, realized
by templatizing the <code class="computeroutput"><span class="keyword">operator</span><span class="special">()</span></code>:
</p>
<p>
</p>
<pre class="programlisting"><span class="keyword">typedef</span> <span class="identifier">boost</span><span class="special">::</span><span class="identifier">numeric</span><span class="special">::</span><span class="identifier">ublas</span><span class="special">::</span><span class="identifier">vector</span><span class="special">&lt;</span> <span class="keyword">double</span> <span class="special">&gt;</span> <span class="identifier">vector_type</span><span class="special">;</span>
<span class="keyword">typedef</span> <span class="identifier">boost</span><span class="special">::</span><span class="identifier">numeric</span><span class="special">::</span><span class="identifier">ublas</span><span class="special">::</span><span class="identifier">matrix</span><span class="special">&lt;</span> <span class="keyword">double</span> <span class="special">&gt;</span> <span class="identifier">matrix_type</span><span class="special">;</span>
<span class="keyword">struct</span> <span class="identifier">stiff_system</span>
<span class="special">{</span>
<span class="keyword">template</span><span class="special">&lt;</span> <span class="keyword">class</span> <span class="identifier">State</span> <span class="special">&gt;</span>
<span class="keyword">void</span> <span class="keyword">operator</span><span class="special">()(</span> <span class="keyword">const</span> <span class="identifier">State</span> <span class="special">&amp;</span><span class="identifier">x</span> <span class="special">,</span> <span class="identifier">State</span> <span class="special">&amp;</span><span class="identifier">dxdt</span> <span class="special">,</span> <span class="keyword">double</span> <span class="identifier">t</span> <span class="special">)</span>
<span class="special">{</span>
<span class="special">...</span>
<span class="special">}</span>
<span class="special">};</span>
<span class="keyword">struct</span> <span class="identifier">stiff_system_jacobi</span>
<span class="special">{</span>
<span class="keyword">template</span><span class="special">&lt;</span> <span class="keyword">class</span> <span class="identifier">State</span> <span class="special">,</span> <span class="keyword">class</span> <span class="identifier">Matrix</span> <span class="special">&gt;</span>
<span class="keyword">void</span> <span class="keyword">operator</span><span class="special">()(</span> <span class="keyword">const</span> <span class="identifier">State</span> <span class="special">&amp;</span><span class="identifier">x</span> <span class="special">,</span> <span class="identifier">Matrix</span> <span class="special">&amp;</span><span class="identifier">J</span> <span class="special">,</span> <span class="keyword">const</span> <span class="keyword">double</span> <span class="special">&amp;</span><span class="identifier">t</span> <span class="special">,</span> <span class="identifier">State</span> <span class="special">&amp;</span><span class="identifier">dfdt</span> <span class="special">)</span>
<span class="special">{</span>
<span class="special">...</span>
<span class="special">}</span>
<span class="special">};</span>
</pre>
<p>
</p>
<p>
Now you can use <code class="computeroutput"><span class="identifier">stiff_system</span></code>
in combination with <code class="computeroutput"><span class="identifier">std</span><span class="special">::</span><span class="identifier">vector</span></code> or <code class="computeroutput"><span class="identifier">boost</span><span class="special">::</span><span class="identifier">array</span></code>.
In the example the explicit time derivative of <span class="emphasis"><em>f(x,t)</em></span>
is introduced separately in the Jacobian. If <span class="emphasis"><em>df / dt = 0</em></span>
simply fill <code class="computeroutput"><span class="identifier">dfdt</span></code> with zeros.
</p>
<p>
A well know solver for stiff systems is the so called Rosenbrock method.
It has a step size control and dense output facilities and can be used like
all the other stepper:
</p>
<p>
</p>
<pre class="programlisting"><span class="identifier">vector_type</span> <span class="identifier">x</span><span class="special">(</span> <span class="number">2</span> <span class="special">,</span> <span class="number">1.0</span> <span class="special">);</span>
<span class="identifier">size_t</span> <span class="identifier">num_of_steps</span> <span class="special">=</span> <span class="identifier">integrate_const</span><span class="special">(</span> <span class="identifier">make_dense_output</span><span class="special">&lt;</span> <span class="identifier">rosenbrock4</span><span class="special">&lt;</span> <span class="keyword">double</span> <span class="special">&gt;</span> <span class="special">&gt;(</span> <span class="number">1.0e-6</span> <span class="special">,</span> <span class="number">1.0e-6</span> <span class="special">)</span> <span class="special">,</span>
<span class="identifier">make_pair</span><span class="special">(</span> <span class="identifier">stiff_system</span><span class="special">()</span> <span class="special">,</span> <span class="identifier">stiff_system_jacobi</span><span class="special">()</span> <span class="special">)</span> <span class="special">,</span>
<span class="identifier">x</span> <span class="special">,</span> <span class="number">0.0</span> <span class="special">,</span> <span class="number">50.0</span> <span class="special">,</span> <span class="number">0.01</span> <span class="special">,</span>
<span class="identifier">cout</span> <span class="special">&lt;&lt;</span> <span class="identifier">phoenix</span><span class="special">::</span><span class="identifier">arg_names</span><span class="special">::</span><span class="identifier">arg2</span> <span class="special">&lt;&lt;</span> <span class="string">" "</span> <span class="special">&lt;&lt;</span> <span class="identifier">phoenix</span><span class="special">::</span><span class="identifier">arg_names</span><span class="special">::</span><span class="identifier">arg1</span><span class="special">[</span><span class="number">0</span><span class="special">]</span> <span class="special">&lt;&lt;</span> <span class="string">"\n"</span> <span class="special">);</span>
</pre>
<p>
</p>
<p>
During the integration 71 steps have been done. Comparing to a classical
Runge-Kutta solver this is a very good result. For example the Dormand-Prince
5 method with step size control and dense output yields 1531 steps.
</p>
<p>
</p>
<pre class="programlisting"><span class="identifier">vector_type</span> <span class="identifier">x2</span><span class="special">(</span> <span class="number">2</span> <span class="special">,</span> <span class="number">1.0</span> <span class="special">);</span>
<span class="identifier">size_t</span> <span class="identifier">num_of_steps2</span> <span class="special">=</span> <span class="identifier">integrate_const</span><span class="special">(</span> <span class="identifier">make_dense_output</span><span class="special">&lt;</span> <span class="identifier">runge_kutta_dopri5</span><span class="special">&lt;</span> <span class="identifier">vector_type</span> <span class="special">&gt;</span> <span class="special">&gt;(</span> <span class="number">1.0e-6</span> <span class="special">,</span> <span class="number">1.0e-6</span> <span class="special">)</span> <span class="special">,</span>
<span class="identifier">stiff_system</span><span class="special">()</span> <span class="special">,</span> <span class="identifier">x2</span> <span class="special">,</span> <span class="number">0.0</span> <span class="special">,</span> <span class="number">50.0</span> <span class="special">,</span> <span class="number">0.01</span> <span class="special">,</span>
<span class="identifier">cout</span> <span class="special">&lt;&lt;</span> <span class="identifier">phoenix</span><span class="special">::</span><span class="identifier">arg_names</span><span class="special">::</span><span class="identifier">arg2</span> <span class="special">&lt;&lt;</span> <span class="string">" "</span> <span class="special">&lt;&lt;</span> <span class="identifier">phoenix</span><span class="special">::</span><span class="identifier">arg_names</span><span class="special">::</span><span class="identifier">arg1</span><span class="special">[</span><span class="number">0</span><span class="special">]</span> <span class="special">&lt;&lt;</span> <span class="string">"\n"</span> <span class="special">);</span>
</pre>
<p>
</p>
<p>
Note, that we have used <a href="http://www.boost.org/doc/libs/1_46_1/libs/spirit/phoenix/doc/html/index.html" target="_top">Boost.Phoenix</a>,
a great functional programming library to create and compose the observer.
</p>
<p>
The full example can be found here: <a href="https://github.com/headmyshoulder/odeint-v2/tree/master/libs/numeric/odeint/examples/stiff_system.cpp" target="_top">stiff_system.cpp</a>
</p>
</div>
<table xmlns:rev="http://www.cs.rpi.edu/~gregod/boost/tools/doc/revision" width="100%"><tr>
<td align="left"></td>
<td align="right"><div class="copyright-footer">Copyright &#169; 2009-2011 Karsten Ahnert
and Mario Mulansky<p>
Distributed under the Boost Software License, Version 1.0. (See accompanying
file LICENSE_1_0.txt or copy at <a href="http://www.boost.org/LICENSE_1_0.txt" target="_top">http://www.boost.org/LICENSE_1_0.txt</a>)
</p>
</div></td>
</tr></table>
<hr>
<div class="spirit-nav">
<a accesskey="p" href="chaotic_systems_and_lyapunov_exponents.html"><img src="../../images/prev.png" alt="Prev"></a><a accesskey="u" href="../tutorial.html"><img src="../../images/up.png" alt="Up"></a><a accesskey="h" href="../../index.html"><img src="../../images/home.png" alt="Home"></a><a accesskey="n" href="special_topics.html"><img src="../../images/next.png" alt="Next"></a>
</div>
</body>
</html>
Reference in New Issue View Git Blame Copy Permalink