Xcos Solvers : Scilab Documentation
[scilab.git] / scilab / modules / xcos / help / en_US / solvers / CVode.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="CVode">
13     <refnamediv>
14         <refname>CVode</refname>
15         <refpurpose>
16             <emphasis>CVode</emphasis> is a numerical solver providing an efficient and stable method to solve ODE Initial Value Problems. Called by scicos(), it uses either <emphasis>BDF</emphasis> or <emphasis>Adams</emphasis> as implicit integration method, and <emphasis>Newton</emphasis> or <emphasis>Functional</emphasis> iterations
17         </refpurpose>
18     </refnamediv>
19     <refsection>
20         <title>Description</title>
21         <para>
22             <emphasis>CVode</emphasis> is a numerical solver providing an efficient and stable 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             Starting with 
33             <emphasis>
34                 y<subscript>0</subscript>
35             </emphasis>
36             ,<emphasis>CVode</emphasis> approximates 
37             <emphasis>
38                 y<subscript>n+1</subscript>
39             </emphasis>
40             with the formula :
41         </para>
42         <para>
43             <latex>
44                 \begin{eqnarray}
45                 \sum_{i=0}^{K_1} \alpha_{n,i} y_{n-i} + h_n\sum_{i=0}^{K_2} \beta_{n,i} \dot{y}_{n-i} = 0,\hspace{10 mm} (1)
46                 \end{eqnarray}
47             </latex>
48             <para>
49                 &#160; with 
50                 <emphasis>
51                     y<subscript>n</subscript>
52                 </emphasis>
53                 the approximation of 
54                 <emphasis>
55                     y(t<subscript>n</subscript>)
56                 </emphasis>
57                 ,and 
58                 <emphasis>
59                     h<subscript>n</subscript>
60                 </emphasis>
61                 =
62                 <emphasis>
63                     t<subscript>n</subscript> - t<subscript>n-1</subscript>
64                 </emphasis>
65                 the step size.
66             </para>
67         </para>
68         <para>
69             These implicit methods are characterized by their respective order <emphasis>q</emphasis>, which indicates the number of intermediate points required to compute 
70             <emphasis>
71                 y<subscript>n+1</subscript>
72             </emphasis>
73             .
74         </para>
75         <para>
76             This is where the difference between <emphasis>BDF</emphasis> and <emphasis>Adams</emphasis> intervenes (<emphasis>Backward Differenciation Formula</emphasis> and <emphasis>Adams-Moulton formula</emphasis>) :
77         </para>
78         <para>
79             <caution>
80                 If the problem is stiff, the user should select <emphasis>BDF</emphasis> :
81             </caution>
82         </para>
83         <itemizedlist>
84             <listitem>
85                 <emphasis>q</emphasis>, the order of the method, is set between 1 and 5 (automated),
86             </listitem>
87             <listitem>
88                 <emphasis>K1 = q</emphasis> and <emphasis>K2 = 0</emphasis>.
89             </listitem>
90         </itemizedlist>
91         <para>
92             In the case of nonstiffness, <emphasis>Adams</emphasis> is preferred :
93         </para>
94         <itemizedlist>
95             <listitem>
96                 <emphasis>q</emphasis> is set between 1 and 12 (automated),
97             </listitem>
98             <listitem>
99                 <emphasis>K1 = 1</emphasis> and <emphasis>K2 = q</emphasis>.
100             </listitem>
101         </itemizedlist>
102         <para>
103             The coefficients are fixed, uniquely determined by the method type, its order, the history of the step sizes, and the normalization 
104             <emphasis>
105                 &#x3B1;<subscript>n, 0</subscript> = -1
106             </emphasis>
107             .
108         </para>
109         <para>
110             For either choice and at each step, injecting this integration in <emphasis>(1)</emphasis> yields the nonlinear system :
111         </para>
112         <para>
113             <latex>
114                 G(y_n)\equiv y_n-h_n\beta_{n,0}f(t_n,y_n)-a_n=0, \hspace{2 mm} where \hspace{2 mm} a_n\equiv \sum_{i>0} (\alpha_{n,i} y_{n-i} + h_n\beta_{n,i}\dot{y}_{n-i})
115             </latex>
116         </para>
117         <para>
118             This system can be solved by either <emphasis>Functional</emphasis> or <emphasis>Newton</emphasis> iterations, described hereafter.
119         </para>
120         <para>
121             In both following cases, the initial "predicted" 
122             <emphasis>
123                 y<subscript>n(0)</subscript>
124             </emphasis>
125             is explicitly computed from the history data, by adding derivatives.
126         </para>
127         <itemizedlist>
128             <para>  </para>
129             <listitem>
130                 <emphasis role="bold">Functional :</emphasis> this method only involves evaluations of <emphasis>f</emphasis>, it simply computes 
131                 <emphasis>
132                     y<subscript>n(0)</subscript>
133                 </emphasis>
134                 by iterating the formula :
135                 <para>
136                     <latex>
137                         y_{n(m+1)} = h_n &#x3B2;_{n,0} f(t_n,y_{n(m+1)}) + a_n
138                     </latex>
139                     <latex>
140                         where \hspace{2 mm} a_n\equiv \sum_{i>0} (\alpha_{n,i} y_{n-i} + h_n\beta_{n,i}\dot{y}_{n-i})
141                     </latex>
142                 </para>
143             </listitem>
144             <para>  </para>
145             <listitem>
146                 <emphasis role="bold">Newton :</emphasis> here, we use an implemented direct dense solver on the linear system :
147                 <para>
148                     <latex>
149                         M[y_{n(m+1)}-y_{n(m)}]=-G(y_{n(m)}), \hspace{4 mm} M \approx I-\gamma J, \hspace{2 mm} J=\frac{\partial f}{\partial y}, \hspace{2 mm} and \hspace{2 mm} \gamma = h_n\beta_{n,0}
150                     </latex>
151                 </para>
152             </listitem>
153         </itemizedlist>
154         <para>
155             In both situations, <emphasis>CVode</emphasis> uses the history array to control the local error 
156             <emphasis>
157                 y<subscript>n(m)</subscript> - y<subscript>n(0)</subscript>
158             </emphasis>
159             and recomputes 
160             <emphasis>
161                 h<subscript>n</subscript>
162             </emphasis>
163             if that error is not satisfying.
164         </para>
165         <para>
166             <emphasis role="bold">
167                 The recommended choices are <emphasis>BDF / Newton</emphasis> for stiff problems and <emphasis>Adams / Functional</emphasis> for the nonstiff ones.
168             </emphasis>
169         </para>
170         <para>
171             The function is called in between activations, because a discrete activation may change the system.
172         </para>
173         <para>
174             Following the criticality of the event (its effect on the continuous problem), we either relaunch the solver with different start and final times as if nothing happened, or, if the system has been modified, we need to "cold-restart" the problem by reinitializing it anew and relaunching the solver.
175         </para>
176     </refsection>
177     <refsection>
178         <title>Examples</title>
179         <para>
180             <link type="scilab" linkend="scilab.xcos/xcos/examples/solvers/ODE_Example.xcos">
181                 <inlinemediaobject>
182                     <imageobject>
183                         <imagedata align="center" fileref="../../../examples/solvers/ODE_Example.xcos" valign="middle"/>
184                     </imageobject>
185                 </inlinemediaobject>
186             </link>
187             <scilab:image><![CDATA[
188 importXcosDiagram(SCI + "/modules/xcos/examples/solvers/ODE_Example.xcos");
189 xcos_simulate(scs_m, 4);
190 ]]></scilab:image>
191         </para>
192         <para>
193             The integral block returns its continuous state, we can evaluate it with BDF / Newton by running the example :
194         </para>
195         <para>
196             <programlisting language="example"><![CDATA[
197       // Import the diagram and set the ending time
198       importXcosDiagram("SCI/modules/xcos/examples/solvers/ODE_Example.xcos");
199       scs_m.props.tf = 3500;
200
201       // Select the solver BDF / Newton
202       scs_m.props.tol(6) = 0;
203
204       // Start the timer, launch the simulation and display time
205       timer();
206       xcos_simulate(scs_m, 4);
207       t = timer();
208       disp(t, "Time for BDF / Newton :");
209       ]]></programlisting>
210         </para>
211         <para>
212             The Scilab console displays :
213             <screen><![CDATA[
214 Time for BDF / Newton :
215  27.549
216             ]]></screen>
217         </para>
218         <para>
219             Now, in the following script, we compare the time difference between the methods by running the example with the four solvers in turn :
220             <link type="scilab" linkend ="scilab.scinotes/xcos/examples/solvers/integ.sce">
221                 Open the script
222             </link>
223         </para>
224         <para>
225             Results :
226             <screen language="example"><![CDATA[
227 Time for BDF / Newton :
228  25.457
229
230 Time for BDF / Functional :
231  24.893
232
233 Time for Adams / Functional :
234  20.049
235
236 Time for Adams / Newton :
237  20.469
238             ]]></screen>
239         </para>
240         <para>
241             The results show that for a simple nonstiff continuous problem, Adams / Functional is fastest.
242         </para>
243     </refsection>
244     <refsection>
245         <title>See Also</title>
246         <simplelist type="inline">
247             <member>
248                 <link linkend="IDA">IDA</link>
249             </member>
250             <member>
251                 <link linkend="RK">Runge-Kutta 4(5)</link>
252             </member>
253             <member>
254                 <link linkend="ode">ode</link>
255             </member>
256             <member>
257                 <link linkend="ode_discrete">ode_discrete</link>
258             </member>
259             <member>
260                 <link linkend="ode_root">ode_root</link>
261             </member>
262             <member>
263                 <link linkend="odedc">odedc</link>
264             </member>
265             <member>
266                 <link linkend="impl">impl</link>
267             </member>
268         </simplelist>
269     </refsection>
270     <refsection>
271         <title>Bibliography</title>
272         <para>
273             <ulink url="https://computation.llnl.gov/casc/sundials/documentation/documentation.html">Sundials Documentation</ulink>
274         </para>
275     </refsection>
276 </refentry>