2cdec40ef293508929aba9207ab278196c68a5c1
[scilab.git] / scilab / modules / xcos / help / en_US / solvers / 3-Dormand-Prince.xml
1 <?xml version="1.0" encoding="UTF-8"?>
2 <!--
3  * Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
4  * Copyright (C) Scilab Enterprises - 2012 - Paul Bignier
5  *
6  * This file must be used under the terms of the CeCILL.
7  * This source file is licensed as described in the file COPYING, which
8  * you should have received as part of this distribution.
9  * The terms are also available at
10  * http://www.cecill.info/licences/Licence_CeCILL_V2.1-en.txt
11  -->
12 <refentry xmlns="http://docbook.org/ns/docbook" xmlns:xlink="http://www.w3.org/1999/xlink" xmlns:svg="http://www.w3.org/2000/svg"  xmlns:mml="http://www.w3.org/1998/Math/MathML" xmlns:db="http://docbook.org/ns/docbook" xmlns:scilab="http://www.scilab.org" xml:lang="en_US" xml:id="DoPri">
13     <refnamediv>
14         <refname>Dormand-Prince 4(5)</refname>
15         <refpurpose>
16             <emphasis>Dormand-Prince</emphasis> is a numerical solver providing an efficient explicit method to solve Ordinary Differential Equations (ODEs) Initial Value Problems.
17         </refpurpose>
18     </refnamediv>
19     <refsection>
20         <title>Description</title>
21         <para>
22             Called by <link linkend="xcos">xcos</link>, <emphasis>Dormand-Prince</emphasis> is a numerical solver providing an efficient fixed-size step method to solve Initial Value Problems of the form:
23         </para>
24         <para>
25             <latex>
26                 \begin{eqnarray}
27                 \dot{y} = f(t,y), \hspace{3 mm} y(t_0) = y_0, \hspace{3 mm} y \in R^N
28                 \end{eqnarray}
29             </latex>
30         </para>
31         <para>
32             <emphasis>CVode</emphasis> and <emphasis>IDA</emphasis> use variable-size steps for the integration.
33         </para>
34         <para>
35             A drawback of that is the unpredictable computation time. With <emphasis>Dormand-Prince</emphasis>, we do not adapt to the complexity of the problem, but we guarantee a stable computation time.
36         </para>
37         <para>
38             As of now, this method is explicit, so it is not concerned with Newton or Functional iterations, and not advised for stiff problems.
39         </para>
40         <para>
41             It is an enhancement of the Euler method, which approximates 
42             <emphasis>
43                 y<subscript>n+1</subscript>
44             </emphasis>
45             by truncating the Taylor expansion.
46         </para>
47         <para>
48             By convention, to use fixed-size steps, the program first computes a fitting <emphasis>h</emphasis> that approaches the simulation parameter <link linkend="Simulatemenu_Menu_entries">max step size</link>.
49         </para>
50         <para>
51             An important difference of <emphasis>Dormand-Prince</emphasis> with the previous methods is that it computes up to the seventh derivative of <emphasis>y</emphasis>, while the others only use linear combinations of <emphasis>y</emphasis> and <emphasis>y'</emphasis>.
52         </para>
53         <para>
54             Here, the next value is determined by the present value 
55             <emphasis>
56                 y<subscript>n</subscript>
57             </emphasis>
58             plus the weighted average of six increments, where each increment is the product of the size of the interval, <emphasis>h</emphasis>, and an estimated slope specified by the function <emphasis>f(t,y)</emphasis>:
59             <itemizedlist>
60                 <listitem>
61                     <emphasis>k1</emphasis> is the increment based on the slope at the beginning of the interval, using 
62                     <emphasis>
63                         y<subscript>n</subscript>
64                     </emphasis>
65                     (Euler's method),
66                 </listitem>
67                 <listitem>
68                     <emphasis>k2, k3, k4</emphasis> and <emphasis>k5</emphasis> are the increments based on the slope at respectively <emphasis>0.2, 0.3, 0.8</emphasis> and <emphasis>0.9</emphasis> of the interval, using combinations of each other,
69                 </listitem>
70                 <listitem>
71                     <emphasis>k6</emphasis> is the increment based on the slope at the end, also using combinations of the other <emphasis>ki</emphasis>.
72                 </listitem>
73             </itemizedlist>
74         </para>
75         <para>
76             We can see that with the <emphasis>ki</emphasis>, we progress in the derivatives of 
77             <emphasis>
78                 y<subscript>n</subscript>
79             </emphasis>
80             . In the computation of the <emphasis>ki</emphasis>, we deliberately use coefficients that yield an error in 
81             <emphasis>
82                 O(h<superscript>5</superscript>)
83             </emphasis>
84             at every step.
85         </para>
86         <para>
87             So the total error is 
88             <emphasis>
89                 number of steps * O(h<superscript>5</superscript>)
90             </emphasis>
91             . And since <emphasis>number of steps = interval size / h</emphasis> by definition, the total error is in 
92             <emphasis>
93                 O(h<superscript>4</superscript>)
94             </emphasis>
95             .
96         </para>
97         <para>
98             That error analysis baptized the method <emphasis>Dormand-Prince 4(5)</emphasis>: 
99             <emphasis>
100                 O(h<superscript>5</superscript>)
101             </emphasis>
102             per step, 
103             <emphasis>
104                 O(h<superscript>4</superscript>)
105             </emphasis>
106             in total.
107         </para>
108         <para>
109             Althought the solver works fine for <link linkend="Simulatemenu_Menu_entries">max step size</link> up to 
110             <emphasis>
111                 10<superscript>-3</superscript>
112             </emphasis>
113             , rounding errors sometimes come into play as it approaches 
114             <emphasis>
115                 4*10<superscript>-4</superscript>
116             </emphasis>
117             . Indeed, the interval splitting cannot be done properly and we get capricious results.
118         </para>
119     </refsection>
120     <refsection>
121         <title>Examples</title>
122         <para>
123             <link type="scilab" linkend="scilab.xcos/xcos/examples/solvers/ODE_Example.zcos">
124                 <inlinemediaobject>
125                     <imageobject>
126                         <imagedata align="center" fileref="../../../examples/solvers/ODE_Example.zcos" valign="middle"/>
127                     </imageobject>
128                 </inlinemediaobject>
129             </link>
130             <scilab:image><![CDATA[
131 loadScicos();
132 loadXcosLibs();
133 importXcosDiagram(SCI + "/modules/xcos/examples/solvers/ODE_Example.zcos");
134 scs_m.props.tol(6) = 5;
135 scs_m.props.tol(7) = 10^-2;
136 try xcos_simulate(scs_m, 4); catch disp(lasterror()); end
137 ]]></scilab:image>
138         </para>
139         <para>
140             The integral block returns its continuous state, we can evaluate it with <emphasis>Dormand-Prince</emphasis> by running the example:
141         </para>
142         <para>
143             <programlisting language="example"><![CDATA[
144       // Import the diagram and set the ending time
145       loadScicos();
146       loadXcosLibs();
147       importXcosDiagram("SCI/modules/xcos/examples/solvers/ODE_Example.zcos");
148       scs_m.props.tf = 5000;
149
150       // Select the Dormand-Prince solver and set the precision
151       scs_m.props.tol(6) = 5;
152       scs_m.props.tol(7) = 10^-2;
153
154       // Start the timer, launch the simulation and display time
155       tic();
156       try xcos_simulate(scs_m, 4); catch disp(lasterror()); end
157       t = toc();
158       disp(t, "Time for Dormand-Prince:");
159       ]]></programlisting>
160         </para>
161         <para>
162             The Scilab console displays:
163             <screen><![CDATA[
164 Time for Dormand-Prince:
165  9.197
166             ]]></screen>
167         </para>
168         <para>
169             Now, in the following script, we compare the time difference between <emphasis>Dormand-Prince</emphasis> and <emphasis>CVode</emphasis> by running the example with the five solvers in turn:
170             <link type="scilab" linkend ="scilab.scinotes/xcos/examples/solvers/integDoPri.sce">
171                 Open the script
172             </link>
173         </para>
174         <para>
175             <screen><![CDATA[
176 Time for BDF / Newton:
177  18.894
178
179 Time for BDF / Functional:
180  18.382
181
182 Time for Adams / Newton:
183  10.368
184
185 Time for Adams / Functional:
186  9.815
187
188 Time for Dormand-Prince:
189  7.842
190             ]]></screen>
191         </para>
192         <para>
193             These results show that on a nonstiff problem, for relatively same precision required and forcing the same step size, <emphasis>Dormand-Prince</emphasis>'s computational overhead (compared to <emphasis>Runge-Kutta</emphasis>) is significant and is close to <emphasis>Adams/Functional</emphasis>. Its error to the solution is althought much smaller than the regular <emphasis>Runge-Kutta 4(5)</emphasis>, for a small overhead in time.
194         </para>
195         <para>
196             Variable-size step ODE solvers are not appropriate for deterministic real-time applications because the computational overhead of taking a time step varies over the course of an application.
197         </para>
198     </refsection>
199     <refsection>
200         <title>See Also</title>
201         <simplelist type="inline">
202             <member>
203                 <link linkend="LSodar">LSodar</link>
204             </member>
205             <member>
206                 <link linkend="CVode">CVode</link>
207             </member>
208             <member>
209                 <link linkend="IDA">IDA</link>
210             </member>
211             <member>
212                 <link linkend="RK">Runge-Kutta 4(5)</link>
213             </member>
214             <member>
215                 <link linkend="ImpRK">Implicit Runge-Kutta 4(5)</link>
216             </member>
217             <member>
218                 <link linkend="DDaskr">DDaskr</link>
219             </member>
220             <member>
221                 <link linkend="Comparisons">Comparisons</link>
222             </member>
223             <member>
224                 <link linkend="ode">ode</link>
225             </member>
226             <member>
227                 <link linkend="ode_discrete">ode_discrete</link>
228             </member>
229             <member>
230                 <link linkend="ode_root">ode_root</link>
231             </member>
232             <member>
233                 <link linkend="odedc">odedc</link>
234             </member>
235             <member>
236                 <link linkend="impl">impl</link>
237             </member>
238         </simplelist>
239     </refsection>
240     <refsection>
241         <title>Bibliography</title>
242         <para>
243             Journal of Computational and Applied Mathematics, Volume 15, Issue 2, 2 June 1986, Pages 203-211 <ulink url="http://dl.acm.org/citation.cfm?id=9958.9963&amp;coll=DL&amp;dl=GUIDE">Dormand-Prince Method</ulink>
244         </para>
245         <para>
246             <ulink url="https://computation.llnl.gov/casc/sundials/documentation/documentation.html">Sundials Documentation</ulink>
247         </para>
248     </refsection>
249     <refsection>
250         <title>History</title>
251         <revhistory>
252             <revision>
253                 <revnumber>5.4.1</revnumber>
254                 <revdescription>Dormand-Prince 4(5) solver added</revdescription>
255             </revision>
256         </revhistory>
257     </refsection>
258 </refentry>