2e62302257f1c3db79f5d45afeefe951748a97f0
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 - 2013 - 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>
15         <refpurpose>
16             <emphasis>DDaskr</emphasis> (short for Double-precision Differential Algebraic equations system Solver with Krylov method and Rootfinding) is a numerical solver providing an efficient and stable method to solve Differential Algebraic Equations systems (DAEs) Initial Value Problems.
17         </refpurpose>
18     </refnamediv>
19     <refsection>
20         <title>Description</title>
21         <para>
22             Called by <link linkend="xcos">xcos</link>, <emphasis>DDaskr</emphasis> (short for Double-precision Differential Algebraic equations system Solver with Krylov method and Rootfinding) 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                 F(t,y,\dot{y}) = 0, \hspace{2 mm} y(t_0)=y_0, \hspace{2 mm} \dot{y}(t_0)=\dot{y}_0, \hspace{3 mm} y, \hspace{1.5 mm} \dot{y}  \hspace{1.5 mm} and \hspace{1.5 mm} F \in R^N \hspace{10 mm} (1)
28                 \end{eqnarray}
29             </latex>
30         </para>
31         <para>
32         </para>
33         Before solving the problem, <emphasis>DDaskr</emphasis> runs an implemented routine to find consistent values for
34         <emphasis>
35             y<subscript>0</subscript>
36         </emphasis>
37         and
38         <emphasis>
39             yPrime<subscript>0</subscript>
40         </emphasis>
41         .
42         <para>
43             Starting then with those
44             <emphasis>
45                 y<subscript>0</subscript>
46             </emphasis>
47             and
48             <emphasis>
49                 yPrime<subscript>0</subscript>
50             </emphasis>
52             <emphasis>
53                 y<subscript>n+1</subscript>
54             </emphasis>
55             with the <emphasis>BDF</emphasis> formula:
56         </para>
57         <para>
58             <latex>
59                 \begin{eqnarray}
60                 \sum_{i=0}^{q} \alpha_{n,i} y_{n-i} = h_n\dot{y}_{n}
61                 \end{eqnarray}
62             </latex>
63             <para>
65                 <emphasis>
66                     y<subscript>n</subscript>
67                 </emphasis>
68                 the approximation of
69                 <emphasis>
70                     y(t<subscript>n</subscript>)
71                 </emphasis>
72                 ,
73                 <emphasis>
74                     h<subscript>n</subscript>
75                 </emphasis>
76                 =
77                 <emphasis>
78                     t<subscript>n</subscript> - t<subscript>n-1</subscript>
79                 </emphasis>
80                 the step size, and the coefficients are fixed, uniquely determined by the method type, its order <emphasis>q</emphasis> ranging from 1 to 5 and the history of the step sizes.
81             </para>
82         </para>
83         <para>
84             Injecting this formula in <emphasis>(1)</emphasis> yields the system:
85         </para>
86         <para>
87             <latex>
88                 G(y_n) \equiv F \left( t_n, \hspace{1.5mm} y_n, \hspace{1.5mm} h_n^{-1}\sum_{i=0}^{q} \alpha_{n,i} y_{n-i} \right) = 0
89             </latex>
90         </para>
91         <para>
92             Its solving is done through a Newton method, but with either <emphasis role="bold">direct</emphasis> or preconditioned GMRes <emphasis role="bold">Krylov</emphasis> iterations.
93         </para>
94         <itemizedlist>
95             <listitem>
96                 <para>
97                     <emphasis role="bold">Direct</emphasis> iteration: start by rewriting the system into:
98                 </para>
99                 <para>
100                     <latex>
101                         J \left[y_{n(m+1)}-y_{n(m)} \right] = -G(y_{n(m)})
102                     </latex>
103                 </para>
104                 <para>
105                     &#160; with <emphasis>J</emphasis> an approximation of the Jacobian:
106                 </para>
107                 <para>
108                     &#160; <latex>
109                         J = \frac{\partial{G}}{\partial{y}} = \frac{\partial{F}}{\partial{y}}+\alpha\frac{\partial{F}}{\partial{\dot{y}}}, \hspace{4 mm} \alpha = \frac{\alpha_{n,0}}{h_n},
110                     </latex>
111                 </para>
112                 <para>
113                     &#160; <emphasis>&#x3B1;</emphasis> changes whenever the step size or the method order varies.
114                 </para>
115                 <para>
116                     Then, an implemented direct dense solver is used and we go on to the next step.
117                 </para>
118             </listitem>
119             <listitem>
120                 <para>
121                     <emphasis role="bold">GMRes Krylov</emphasis> method: first, precondition the system by applying the Jacobian matrix mentioned above.
122                 </para>
123                 <para>
124                     Secondly, compute the next Krylov space basis and update the Hessenberg matrix.
125                 </para>
126                 <para>
127                     Test for convergence for the first time. If it doesn't pass, calculate the residual, which will lead to a new potential solution, and iterate until that residual satisfies convergence.
128                 </para>
129             </listitem>
130         </itemizedlist>
131         <para>
132             <emphasis>DDaskr</emphasis> uses the history array to control the local error
133             <emphasis>
134                 y<subscript>n(m)</subscript> - y<subscript>n(0)</subscript>
135             </emphasis>
136             and recomputes
137             <emphasis>
138                 h<subscript>n</subscript>
139             </emphasis>
140             if that error is not satisfying.
141         </para>
142         <para>
143             The function is called in between activations, because a discrete activation may change the system.
144         </para>
145         <para>
146             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.
147         </para>
148         <para>
149             Averagely, <emphasis>DDaskr</emphasis> accepts tolerances up to 10<superscript>-11</superscript>. Beyond that, it returns a <emphasis>Too much accuracy</emphasis> requested error.
150         </para>
151         <para>
152             As of now, <emphasis>DDaskr</emphasis> can only be applied on systems where the root functions are punctually 0, and not 0-flat.
153         </para>
154     </refsection>
156         <title>Example</title>
157         <para>
158             The 'Modelica Generic' block returns its continuous states, we can evaluate them with <emphasis>DDaskr</emphasis> by running the example:
159         </para>
160         <para>
162                 <inlinemediaobject>
163                     <imageobject>
164                         <imagedata align="center" fileref="../../../examples/solvers/DAE_Example.zcos" valign="middle"/>
165                     </imageobject>
166                 </inlinemediaobject>
168             <scilab:image localized="true"><![CDATA[
171 importXcosDiagram(SCI + "/modules/xcos/examples/solvers/DAE_Example.zcos");
172 // Redefining messagebox() to avoid popup
173 function messagebox(msg, Title)
174  disp(Title);
175  disp(msg);
176 endfunction
177 scs_m.props.tol(6) = 101;
178 try xcos_simulate(scs_m, 4); catch disp(lasterror()); end
179 ]]></scilab:image>
180         </para>
181     </refsection>
182     <refsection>
183         <title>See Also</title>
184         <simplelist type="inline">
185             <member>
187             </member>
188             <member>
190             </member>
191             <member>
193             </member>
194             <member>
196             </member>
197             <member>
199             </member>
200             <member>
202             </member>
203             <member>
205             </member>
206             <member>
208             </member>
209             <member>
211             </member>
212             <member>
214             </member>
215             <member>
217             </member>
218             <member>
220             </member>
221             <member>
223             </member>
224         </simplelist>
225     </refsection>
226     <refsection>
227         <title>Bibliography</title>
228         <para>
230         </para>
231         <para>