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  -->
13     <refnamediv>
14         <refname>CVode</refname>
15         <refpurpose>
16             <emphasis>CVode</emphasis> (short for C-language Variable-coefficients ODE solver) is a numerical solver providing an efficient and stable method to solve Ordinary Differential Equations (ODEs) Initial Value Problems. 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             Called by <link linkend="xcos">xcos</link>, <emphasis>CVode</emphasis> (short for C-language Variable-coefficients ODE solver) 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                  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 Differentiation 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                 α<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 scilab:localized="true">
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 β_{n,0} f(t_n,y_{n(m)}) + a_n
138                     </latex>
139                     <latex scilab:localized="true">
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 scilab:localized="true">
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         <para>
177             Averagely, <emphasis>CVode</emphasis> accepts tolerances up to 10<superscript>-16</superscript>. Beyond that, it returns a <emphasis>Too much accuracy</emphasis> requested error.
178         </para>
179     </refsection>
180     <refsection>
181         <title>Examples</title>
182         <para>
184                 <inlinemediaobject>
185                     <imageobject>
186                         <imagedata align="center" fileref="../../../examples/solvers/ODE_Example.zcos" valign="middle"/>
187                     </imageobject>
188                 </inlinemediaobject>
190             <scilab:image><![CDATA[
193 importXcosDiagram(SCI + "/modules/xcos/examples/solvers/ODE_Example.zcos");
194 scs_m.props.tol(6) = 1;
195 try xcos_simulate(scs_m, 4); catch disp(lasterror()); end
196 ]]></scilab:image>
197         </para>
198         <para>
199             The integral block returns its continuous state, we can evaluate it with <emphasis>BDF / Newton</emphasis> by running the example:
200         </para>
201         <para>
202             <programlisting language="example"><![CDATA[
203       // Import the diagram and set the ending time
206       importXcosDiagram("SCI/modules/xcos/examples/solvers/ODE_Example.zcos");
207       scs_m.props.tf = 5000;
209       // Select the solver BDF / Newton
210       scs_m.props.tol(6) = 1;
212       // Start the timer, launch the simulation and display time
213       tic();
214       try xcos_simulate(scs_m, 4); catch disp(lasterror()); end
215       t = toc();
216       disp(t, "Time for BDF / Newton:");
217       ]]></programlisting>
218         </para>
219         <para>
220             The Scilab console displays:
221             <screen><![CDATA[
222 Time for BDF / Newton:
223  13.438
224             ]]></screen>
225         </para>
226         <para>
227             Now, in the following script, we compare the time difference between the methods by running the example with the four solvers in turn:
229                 Open the script
231         </para>
232         <para>
233             Results:
234             <screen language="example"><![CDATA[
235 Time for BDF / Newton:
236  18.894
238 Time for BDF / Functional:
239  18.382
241 Time for Adams / Newton:
242  10.368
244 Time for Adams / Functional:
245  9.815
246             ]]></screen>
247         </para>
248         <para>
249             The results show that for a simple nonstiff continuous problem, <emphasis>Adams / Functional</emphasis> is fastest.
250         </para>
251     </refsection>
252     <refsection>
253         <title>See Also</title>
254         <simplelist type="inline">
255             <member>
257             </member>
258             <member>
260             </member>
261             <member>
263             </member>
264             <member>
266             </member>
267             <member>
269             </member>
270             <member>
272             </member>
273             <member>
275             </member>
276             <member>
278             </member>
279             <member>
281             </member>
282             <member>
284             </member>
285             <member>