1<html> 2<head> 3<meta http-equiv="Content-Type" content="text/html; charset=UTF-8"> 4<title>Stiff systems</title> 5<link rel="stylesheet" href="../../../../../../../doc/src/boostbook.css" type="text/css"> 6<meta name="generator" content="DocBook XSL Stylesheets V1.79.1"> 7<link rel="home" href="../../index.html" title="Chapter 1. Boost.Numeric.Odeint"> 8<link rel="up" href="../tutorial.html" title="Tutorial"> 9<link rel="prev" href="chaotic_systems_and_lyapunov_exponents.html" title="Chaotic systems and Lyapunov exponents"> 10<link rel="next" href="complex_state_types.html" title="Complex state types"> 11</head> 12<body bgcolor="white" text="black" link="#0000FF" vlink="#840084" alink="#0000FF"> 13<table cellpadding="2" width="100%"><tr> 14<td valign="top"><img alt="Boost C++ Libraries" width="277" height="86" src="../../logo.jpg"></td> 15<td align="center"><a href="../../../../../../../index.html">Home</a></td> 16<td align="center"><a href="../../../../../../../libs/libraries.htm">Libraries</a></td> 17<td align="center"><a href="http://www.boost.org/users/people.html">People</a></td> 18<td align="center"><a href="http://www.boost.org/users/faq.html">FAQ</a></td> 19<td align="center"><a href="../../../../../../../more/index.htm">More</a></td> 20</tr></table> 21<hr> 22<div class="spirit-nav"> 23<a accesskey="p" href="chaotic_systems_and_lyapunov_exponents.html"><img src="../../../../../../../doc/src/images/prev.png" alt="Prev"></a><a accesskey="u" href="../tutorial.html"><img src="../../../../../../../doc/src/images/up.png" alt="Up"></a><a accesskey="h" href="../../index.html"><img src="../../../../../../../doc/src/images/home.png" alt="Home"></a><a accesskey="n" href="complex_state_types.html"><img src="../../../../../../../doc/src/images/next.png" alt="Next"></a> 24</div> 25<div class="section"> 26<div class="titlepage"><div><div><h3 class="title"> 27<a name="boost_numeric_odeint.tutorial.stiff_systems"></a><a class="link" href="stiff_systems.html" title="Stiff systems">Stiff systems</a> 28</h3></div></div></div> 29<p> 30 An important class of ordinary differential equations are so called stiff 31 system which are characterized by two or more time scales of different order. 32 Examples of such systems are found in chemical systems where reaction rates 33 of individual sub-reaction might differ over large ranges, for example: 34 </p> 35<p> 36 <span class="emphasis"><em>d S<sub>1</sub> / dt = - 101 S<sub>2</sub> - 100 S<sub>1</sub></em></span> 37 </p> 38<p> 39 <span class="emphasis"><em>d S<sub>2</sub> / dt = S<sub>1</sub></em></span> 40 </p> 41<p> 42 In order to efficiently solve stiff systems numerically the Jacobian 43 </p> 44<p> 45 <span class="emphasis"><em>J = d f<sub>i</sub> / d x<sub>j</sub></em></span> 46 </p> 47<p> 48 is needed. Here is the definition of the above example 49 </p> 50<p> 51</p> 52<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"><</span> <span class="keyword">double</span> <span class="special">></span> <span class="identifier">vector_type</span><span class="special">;</span> 53<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"><</span> <span class="keyword">double</span> <span class="special">></span> <span class="identifier">matrix_type</span><span class="special">;</span> 54 55<span class="keyword">struct</span> <span class="identifier">stiff_system</span> 56<span class="special">{</span> 57 <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">&</span><span class="identifier">x</span> <span class="special">,</span> <span class="identifier">vector_type</span> <span class="special">&</span><span class="identifier">dxdt</span> <span class="special">,</span> <span class="keyword">double</span> <span class="comment">/* t */</span> <span class="special">)</span> 58 <span class="special">{</span> 59 <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> 60 <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> 61 <span class="special">}</span> 62<span class="special">};</span> 63 64<span class="keyword">struct</span> <span class="identifier">stiff_system_jacobi</span> 65<span class="special">{</span> 66 <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">&</span> <span class="comment">/* x */</span> <span class="special">,</span> <span class="identifier">matrix_type</span> <span class="special">&</span><span class="identifier">J</span> <span class="special">,</span> <span class="keyword">const</span> <span class="keyword">double</span> <span class="special">&</span> <span class="comment">/* t */</span> <span class="special">,</span> <span class="identifier">vector_type</span> <span class="special">&</span><span class="identifier">dfdt</span> <span class="special">)</span> 67 <span class="special">{</span> 68 <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> 69 <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> 70 <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> 71 <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> 72 <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> 73 <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> 74 <span class="special">}</span> 75<span class="special">};</span> 76</pre> 77<p> 78 </p> 79<p> 80 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> 81 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> 82 since the stiff integrator only accepts these types. However, you might want 83 use non-stiff integrators on this system, too - we will do so later for demonstration. 84 Therefore we want to use the same function also with other state_types, realized 85 by templatizing the <code class="computeroutput"><span class="keyword">operator</span><span class="special">()</span></code>: 86 </p> 87<p> 88</p> 89<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"><</span> <span class="keyword">double</span> <span class="special">></span> <span class="identifier">vector_type</span><span class="special">;</span> 90<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"><</span> <span class="keyword">double</span> <span class="special">></span> <span class="identifier">matrix_type</span><span class="special">;</span> 91 92<span class="keyword">struct</span> <span class="identifier">stiff_system</span> 93<span class="special">{</span> 94 <span class="keyword">template</span><span class="special"><</span> <span class="keyword">class</span> <span class="identifier">State</span> <span class="special">></span> 95 <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">&</span><span class="identifier">x</span> <span class="special">,</span> <span class="identifier">State</span> <span class="special">&</span><span class="identifier">dxdt</span> <span class="special">,</span> <span class="keyword">double</span> <span class="identifier">t</span> <span class="special">)</span> 96 <span class="special">{</span> 97 <span class="special">...</span> 98 <span class="special">}</span> 99<span class="special">};</span> 100 101<span class="keyword">struct</span> <span class="identifier">stiff_system_jacobi</span> 102<span class="special">{</span> 103 <span class="keyword">template</span><span class="special"><</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">></span> 104 <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">&</span><span class="identifier">x</span> <span class="special">,</span> <span class="identifier">Matrix</span> <span class="special">&</span><span class="identifier">J</span> <span class="special">,</span> <span class="keyword">const</span> <span class="keyword">double</span> <span class="special">&</span><span class="identifier">t</span> <span class="special">,</span> <span class="identifier">State</span> <span class="special">&</span><span class="identifier">dfdt</span> <span class="special">)</span> 105 <span class="special">{</span> 106 <span class="special">...</span> 107 <span class="special">}</span> 108<span class="special">};</span> 109</pre> 110<p> 111 </p> 112<p> 113 Now you can use <code class="computeroutput"><span class="identifier">stiff_system</span></code> 114 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>. 115 In the example the explicit time derivative of <span class="emphasis"><em>f(x,t)</em></span> 116 is introduced separately in the Jacobian. If <span class="emphasis"><em>df / dt = 0</em></span> 117 simply fill <code class="computeroutput"><span class="identifier">dfdt</span></code> with zeros. 118 </p> 119<p> 120 A well know solver for stiff systems is the Rosenbrock method. It has a step 121 size control and dense output facilities and can be used like all the other 122 steppers: 123 </p> 124<p> 125</p> 126<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> 127 128<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"><</span> <span class="identifier">rosenbrock4</span><span class="special"><</span> <span class="keyword">double</span> <span class="special">></span> <span class="special">>(</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> 129 <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> 130 <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> 131 <span class="identifier">cout</span> <span class="special"><<</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"><<</span> <span class="string">" "</span> <span class="special"><<</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"><<</span> <span class="string">"\n"</span> <span class="special">);</span> 132</pre> 133<p> 134 </p> 135<p> 136 During the integration 71 steps have been done. Comparing to a classical 137 Runge-Kutta solver this is a very good result. For example the Dormand-Prince 138 5 method with step size control and dense output yields 1531 steps. 139 </p> 140<p> 141</p> 142<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> 143 144<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"><</span> <span class="identifier">runge_kutta_dopri5</span><span class="special"><</span> <span class="identifier">vector_type</span> <span class="special">></span> <span class="special">>(</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> 145 <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> 146 <span class="identifier">cout</span> <span class="special"><<</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"><<</span> <span class="string">" "</span> <span class="special"><<</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"><<</span> <span class="string">"\n"</span> <span class="special">);</span> 147</pre> 148<p> 149 </p> 150<p> 151 Note, that we have used <a href="http://www.boost.org/doc/libs/release/libs/phoenix/" target="_top">Boost.Phoenix</a>, 152 a great functional programming library, to create and compose the observer. 153 </p> 154<p> 155 The full example can be found here: <a href="https://github.com/headmyshoulder/odeint-v2/blob/master/examples/stiff_system.cpp" target="_top">stiff_system.cpp</a> 156 </p> 157</div> 158<table xmlns:rev="http://www.cs.rpi.edu/~gregod/boost/tools/doc/revision" width="100%"><tr> 159<td align="left"></td> 160<td align="right"><div class="copyright-footer">Copyright © 2009-2015 Karsten Ahnert and Mario Mulansky<p> 161 Distributed under the Boost Software License, Version 1.0. (See accompanying 162 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>) 163 </p> 164</div></td> 165</tr></table> 166<hr> 167<div class="spirit-nav"> 168<a accesskey="p" href="chaotic_systems_and_lyapunov_exponents.html"><img src="../../../../../../../doc/src/images/prev.png" alt="Prev"></a><a accesskey="u" href="../tutorial.html"><img src="../../../../../../../doc/src/images/up.png" alt="Up"></a><a accesskey="h" href="../../index.html"><img src="../../../../../../../doc/src/images/home.png" alt="Home"></a><a accesskey="n" href="complex_state_types.html"><img src="../../../../../../../doc/src/images/next.png" alt="Next"></a> 169</div> 170</body> 171</html> 172
