Add sparse unitary test about matrix insertion
[scilab.git] / scilab / modules / xcos / help / en_US / solvers / Dormand-Price.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-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 DoPri, 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 DoPri 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>Dopri 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.xcos">
124                 <inlinemediaobject>
125                     <imageobject>
126                         <imagedata align="center" fileref="../../../examples/solvers/ODE_Example.xcos" 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.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 DoPri 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.xcos");
148       scs_m.props.tf = 10000;
149
150       // Select the DoPri 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 DoPri :");
159       ]]></programlisting>
160         </para>
161         <para>
162             The Scilab console displays :
163             <screen><![CDATA[
164 Time for DoPri :
165  31.501
166             ]]></screen>
167         </para>
168         <para>
169             Now, in the following script, we compare the time difference between DoPri and Sundials 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  37.882
178
179 Time for BDF / Functional :
180  32.815
181
182 Time for Adams / Newton :
183  30.108
184
185 Time for Adams / Functional :
186  29.242
187
188 Time for DoPri :
189  31.501
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, DoPri's computational overhead is significant. 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>
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="ode">ode</link>
216             </member>
217             <member>
218                 <link linkend="ode_discrete">ode_discrete</link>
219             </member>
220             <member>
221                 <link linkend="ode_root">ode_root</link>
222             </member>
223             <member>
224                 <link linkend="odedc">odedc</link>
225             </member>
226             <member>
227                 <link linkend="impl">impl</link>
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>
237             <ulink url="https://computation.llnl.gov/casc/sundials/documentation/documentation.html">Sundials Documentation</ulink>
238         </para>
239     </refsection>
240     <refsection>
241         <title>History</title>
242         <revhistory>
243             <revision>
244                 <revnumber>5.4.1</revnumber>
245                 <revdescription>Dormand-Price 4(5) solver added</revdescription>
246             </revision>
247         </revhistory>
248     </refsection>
249 </refentry>