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-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-Price 4(5)</refname>
15         <refpurpose>
16             <emphasis>Dormand-Price</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-Price</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 Dormand-Price, 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 Dormand-Price 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-Price 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>
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>
124                 <inlinemediaobject>
125                     <imageobject>
126                         <imagedata align="center" fileref="../../../examples/solvers/ODE_Example.xcos" valign="middle"/>
127                     </imageobject>
128                 </inlinemediaobject>
130             <scilab:image><![CDATA[
133 importXcosDiagram(SCI + "/modules/xcos/examples/solvers/ODE_Example.xcos");
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 Dormand-Price by running the example :
141         </para>
142         <para>
143             <programlisting language="example"><![CDATA[
144       // Import the diagram and set the ending time
147       importXcosDiagram("SCI/modules/xcos/examples/solvers/ODE_Example.xcos");
148       scs_m.props.tf = 5000;
150       // Select the Dormand-Price solver and set the precision
151       scs_m.props.tol(6) = 5;
152       scs_m.props.tol(7) = 10^-2;
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-Price :");
159       ]]></programlisting>
160         </para>
161         <para>
162             The Scilab console displays :
163             <screen><![CDATA[
164 Time for Dormand-Price :
165  17.126
166             ]]></screen>
167         </para>
168         <para>
169             Now, in the following script, we compare the time difference between Dormand-Price and Sundials by running the example with the five solvers in turn :
171                 Open the script
173         </para>
174         <para>
175             <screen><![CDATA[
176 Time for BDF / Newton :
177  12.852
179 Time for BDF / Functional :
180  12.343
182 Time for Adams / Newton :
183  8.75
185 Time for Adams / Functional :
186  7.789
188 Time for DoPri :
189  7.606
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, Dormand-Price's computational overhead (compared to Runge-Kutta) is significant and is close to <emphasis>Adams/Functional</emphasis>. Its error to the solution is althought much smaller than the regular Runge-Kutta 4(5), for a small overhead in time.
194         </para>
195         <para>
196             Variable step-size 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>
204             </member>
205             <member>
207             </member>
208             <member>
210             </member>
211             <member>
213             </member>
214             <member>
216             </member>
217             <member>
219             </member>
220             <member>
222             </member>
223             <member>
225             </member>
226             <member>
228             </member>
229         </simplelist>
230     </refsection>
231     <refsection>
232         <title>Bibliography</title>
233         <para>
234             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-Price Method</ulink>
235         </para>
236         <para>