Differential_equations doc: warning fix
[scilab.git] / scilab / modules / differential_equations / help / en_US / ode.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) 2008 - INRIA
5  * Copyright (C) 2011 - DIGITEO - Michael Baudin
6  * Copyright (C) 2012 - Scilab Enterprises - Adeline CARNIS
7  *
8  * This file must be used under the terms of the CeCILL.
9  * This source file is licensed as described in the file COPYING, which
10  * you should have received as part of this distribution.  The terms
11  * are also available at
12  * http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt
13  *
14  -->
15 <refentry xmlns="http://docbook.org/ns/docbook" xmlns:xlink="http://www.w3.org/1999/xlink" xmlns:svg="http://www.w3.org/2000/svg" xmlns:ns5="http://www.w3.org/1999/xhtml" xmlns:mml="http://www.w3.org/1998/Math/MathML" xmlns:db="http://docbook.org/ns/docbook" xmlns:scilab="http://www.scilab.org" xml:id="ode" xml:lang="en">
16     <refnamediv>
17         <refname>ode</refname>
18         <refpurpose>ordinary differential equation solver</refpurpose>
19     </refnamediv>
20     <refsynopsisdiv>
21         <title>Calling Sequence</title>
22         <synopsis>
23             y = ode(y0, t0, t, f)
24             [y, w, iw] = ode([type,] y0, t0, t [,rtol [,atol]], f [,jac] [,w, iw])
25             [y, rd, w, iw] = ode("root", y0, t0, t [,rtol [,atol]], f [,jac],ng, g [,w,iw])
26             y = ode("discrete", y0, k0, kvect, f)
27         </synopsis>
28     </refsynopsisdiv>
29     <refsection>
30         <title>Arguments</title>
31         <variablelist>
32             <varlistentry>
33                 <term>y0</term>
34                 <listitem>
35                     <para>a real vector or matrix, the initial conditions.</para>
36                 </listitem>
37             </varlistentry>
38             <varlistentry>
39                 <term>t0</term>
40                 <listitem>
41                     <para>
42                         a real scalar, the initial time.
43                     </para>
44                 </listitem>
45             </varlistentry>
46             <varlistentry>
47                 <term>t</term>
48                 <listitem>
49                     <para>
50                         a real vector, the times at which the solution is computed.
51                     </para>
52                 </listitem>
53             </varlistentry>
54             <varlistentry>
55                 <term>f</term>
56                 <listitem>
57                     <para>
58                         a function, external, string or list, the right hand side of
59                         the differential equation.
60                     </para>
61                 </listitem>
62             </varlistentry>
63             <varlistentry>
64                 <term>type</term>
65                 <listitem>
66                     <para>
67                         a string, the solver to use. The available solvers are
68                         <literal>"adams"</literal>,
69                         <literal>"stiff"</literal>, <literal>"rk"</literal>,
70                         <literal>"rkf"</literal>, <literal>"fix"</literal>,
71                         <literal>"discrete"</literal>, <literal>"roots"</literal>.
72                     </para>
73                 </listitem>
74             </varlistentry>
75             <varlistentry>
76                 <term>rtol</term>
77                 <listitem>
78                     <para>
79                         a real constant or real vector of the same size as
80                         <varname>y</varname>, the relative tolerance.
81                     </para>
82                 </listitem>
83             </varlistentry>
84             <varlistentry>
85                 <term>atol</term>
86                 <listitem>
87                     <para>
88                         a real constant or real vectors of the same size as
89                         <varname>y</varname>, the absolute tolerance.
90                     </para>
91                 </listitem>
92             </varlistentry>
93             <varlistentry>
94                 <term>jac</term>
95                 <listitem>
96                     <para>
97                         a function, external, string or list, the Jacobian of the
98                         function <varname>f</varname>.
99                     </para>
100                 </listitem>
101             </varlistentry>
102             <varlistentry>
103                 <term>w, iw</term>
104                 <listitem>
105                     <para>real vectors. (INPUT/OUTPUT)</para>
106                 </listitem>
107             </varlistentry>
108             <varlistentry>
109                 <term>ng</term>
110                 <listitem>
111                     <para>an integer.</para>
112                 </listitem>
113             </varlistentry>
114             <varlistentry>
115                 <term>g</term>
116                 <listitem>
117                     <para>an external (function or character string or list).</para>
118                 </listitem>
119             </varlistentry>
120             <varlistentry>
121                 <term>k0</term>
122                 <listitem>
123                     <para>an integer (initial time).</para>
124                 </listitem>
125             </varlistentry>
126             <varlistentry>
127                 <term>kvect</term>
128                 <listitem>
129                     <para>an integer vector.</para>
130                 </listitem>
131             </varlistentry>
132             <varlistentry>
133                 <term>y</term>
134                 <listitem>
135                     <para>a real vector or matrix. (OUTPUT)</para>
136                 </listitem>
137             </varlistentry>
138             <varlistentry>
139                 <term>rd</term>
140                 <listitem>
141                     <para>a real vector. (OUTPUT)</para>
142                 </listitem>
143             </varlistentry>
144         </variablelist>
145     </refsection>
146     <refsection>
147         <title>Description</title>
148         <para>
149             <function>ode</function> solves explicit Ordinary Different Equations defined by:
150         </para>
151         <para>
152             <latex>
153                 \begin{eqnarray}
154                 \begin{array}{l}
155                 \dfrac{dy}{dt} = f(t,y)\\
156                 y(t_0)=y_0
157                 \end{array}
158                 \end{eqnarray}
159             </latex>
160         </para>
161         <para>
162             It is an interface to various solvers, in particular to ODEPACK.
163         </para>
164         <para>
165             In this help, we only describe the use of <function>ode</function> for
166             standard explicit ODE systems.
167         </para>
168         <para>
169             The simplest call of <function>ode</function> is:
170             <literal>y = ode(y0,t0,t,f)</literal> where <literal>y0</literal> is the
171             vector of initial conditions, <literal>t0</literal> is the initial
172             time, <literal>t</literal> is the vector of times at which the
173             solution <varname>y</varname> is computed and <varname>y</varname> is
174             matrix of solution vectors
175             <literal>y=[y(t(1)),y(t(2)),...]</literal>.
176         </para>
177         <para>
178             The input argument <varname>f</varname> defines the right hand side of the
179             first order differential equation. This argument is a function with a specific header.
180         </para>
181         <itemizedlist>
182             <listitem>
183                 <para>
184                     If <varname>f</varname> is a Scilab function, its calling sequence
185                     must be
186                 </para>
187                 <screen><![CDATA[
188 ydot = f(t,y)
189  ]]></screen>
190                 <para>
191                     where <literal>t</literal> is a real scalar (the time) and
192                     <varname>y</varname> is a real vector (the state) and
193                     <literal>ydot</literal> is a real vector (the first order derivative
194                     dy/dt).
195                 </para>
196             </listitem>
197             <listitem>
198                 <para>
199                     If <varname>f</varname> is a string, it is
200                     the name of a Fortran subroutine or a C compiled function. For example, if
201                     we call <literal>ode(y0,t0,t,"fex")</literal>, then the
202                     subroutine <literal>fex</literal> is called.
203                 </para>
204                 <para>
205                     The Fortran routine must have the header:
206                 </para>
207                 <screen><![CDATA[
208 fex(n,t,y,ydot)
209  ]]></screen>
210                 <para>
211                     where <literal>n</literal> is an integer, <literal>t</literal> is
212                     a double precision scalar, <varname>y</varname> and
213                     <literal>ydot</literal> are double precision
214                     vectors.
215                 </para>
216                 <para>
217                     The C function must have the header:
218                 </para>
219                 <screen><![CDATA[
220 fex(int *n,double *t,double *y,double *ydot)
221  ]]></screen>
222                 <para>
223                     where <literal>t</literal> is the time, <varname>y</varname> the
224                     state and <literal>ydot</literal> is the state derivative
225                     (dy/dt).
226                 </para>
227                 <para>
228                     This external can be build in a OS independent way using
229                     <link linkend="ilib_for_link">ilib_for_link</link> and dynamically
230                     linked to Scilab by the <link linkend="link">link</link>
231                     function.
232                 </para>
233             </listitem>
234             <listitem>
235                 <para>
236                     It may happen that the simulator <varname>f</varname> needs extra
237                     arguments.
238                     In this case, we can use the following feature.
239                     The <varname>f</varname> argument can also be a list
240                     <literal>lst=list(simuf,u1,u2,...un)</literal> where
241                     <literal>simuf</literal> is a Scilab function with syntax:
242                     <literal>ydot = f(t,y,u1,u2,...,un)</literal>
243                     and <literal>u1</literal>, <literal>u2</literal>, ...,
244                     <literal>un</literal> are extra arguments which are automatically
245                     passed to the simulator <literal>simuf</literal>.
246                 </para>
247             </listitem>
248         </itemizedlist>
249         <para>
250             The function <varname>f</varname> can return a p-by-q matrix instead of a vector.
251             With this matrix notation, we
252             solve the <literal>n=p+q</literal> ODE's system
253             <literal>dY/dt=F(t,Y)</literal> where <literal>Y</literal> is a
254             <literal>p x q</literal> matrix. Then initial conditions,
255             <literal>Y0</literal>, must also be a <literal>p x q</literal> matrix
256             and the result of <function>ode</function> is the <literal>p-by-q(T+1)</literal> matrix
257             <literal>[Y(t_0),Y(t_1),...,Y(t_T)]</literal>.
258         </para>
259         <para>
260             The tolerances <varname>rtol</varname> and <varname>atol</varname> are
261             thresholds for relative and absolute estimated errors. The estimated
262             error on <literal>y(i)</literal> is:
263             <literal>rtol(i)*abs(y(i))+atol(i)</literal>
264             and integration is carried out as far as this error is small for
265             all components of the state. If <varname>rtol</varname> and/or
266             <varname>atol</varname> is a constant <literal>rtol(i)</literal>
267             and/or <literal>atol(i)</literal> are set to this constant value.
268             Default values for <varname>rtol</varname> and <varname>atol</varname>
269             are respectively <literal>rtol=1.d-5</literal> and
270             <literal>atol=1.d-7</literal> for most solvers and
271             <literal>rtol=1.d-3</literal> and <literal>atol=1.d-4</literal> for
272             <literal>"rfk"</literal> and <literal>"fix"</literal>.
273         </para>
274         <para>
275             For stiff problems, it is better to give the Jacobian of the RHS
276             function as the optional argument <varname>jac</varname>.
277             The Jacobian  is an external i.e. a function with specified syntax, or the name of a
278             Fortran subroutine or a C function (character string) with specified
279             calling sequence or a list.
280         </para>
281         <itemizedlist>
282             <listitem>
283                 <para>
284                     If <varname>jac</varname> is a function the syntax should be
285                     <literal>J=jac(t,y)</literal> where <literal>t</literal> is a real
286                     scalar (time) and <varname>y</varname> a real vector (state).
287                     The result matrix <literal>J</literal> must evaluate to df/dx i.e.
288                     <literal>J(k,i) = dfk/dxi</literal> where <literal>fk</literal> is the
289                     <literal>k</literal>-th component of <varname>f</varname>.
290                 </para>
291             </listitem>
292             <listitem>
293                 <para>
294                     If <varname>jac</varname> is a character string it refers to the
295                     name of a Fortran subroutine or a C function.
296                 </para>
297                 <para>
298                     The Fortran routine must have the header:
299                 </para>
300                 <screen><![CDATA[
301 subroutine fex(n,t,y,ml,mu,J,nrpd)
302 integer n,ml,mu,nrpd
303 double precision t,y(*),J(*)
304  ]]></screen>
305                 <para>
306                     The C function must have the header:
307                 </para>
308                 <screen><![CDATA[
309 void fex(int *n,double *t,double *y,int *ml,int *mu,double *J,int *nrpd,)
310  ]]></screen>
311                 <para>
312                     In most cases you have not to refer <literal>ml</literal>, <literal>mu</literal> and
313                     <literal>nrpd</literal>.
314                 </para>
315             </listitem>
316             <listitem>
317                 <para>
318                     If <varname>jac</varname> is a list the same conventions as for
319                     <varname>f</varname> apply.
320                 </para>
321             </listitem>
322         </itemizedlist>
323         <para>
324             Optional arguments <varname>w</varname> and
325             <varname>iw</varname> are vectors for storing information returned by
326             the integration routine (see <link linkend="ode_optional_output">ode_optional_output</link> for details).
327             When these vectors are provided in RHS of <function>ode</function> the
328             integration re-starts with the same parameters as in its previous
329             stop.
330         </para>
331         <para>
332             More options can be given to ODEPACK solvers by using
333             <literal>%ODEOPTIONS</literal> variable. See <link linkend="odeoptions">odeoptions</link>.
334         </para>
335     </refsection>
336     <refsection>
337         <title>The solvers</title>
338         <para>
339             The type of problem solved and
340             the method used depend on the value of the first optional argument
341             <varname>type</varname> which can be one of the following strings:
342         </para>
343         <variablelist>
344             <varlistentry>
345                 <term>&lt;not given&gt;:</term>
346                 <listitem>
347                     <para>
348                         <literal>lsoda</literal> solver of package ODEPACK is called
349                         by default. It automatically selects between nonstiff
350                         predictor-corrector Adams method and stiff Backward Differentiation
351                         Formula (BDF) method. It uses nonstiff method initially and
352                         dynamically monitors data in order to decide which method to
353                         use.
354                     </para>
355                 </listitem>
356             </varlistentry>
357             <varlistentry>
358                 <term>"adams":</term>
359                 <listitem>
360                     <para>
361                         This is for nonstiff problems. <literal>lsode</literal> solver
362                         of package ODEPACK is called and it uses the Adams method.
363                     </para>
364                 </listitem>
365             </varlistentry>
366             <varlistentry>
367                 <term>"stiff":</term>
368                 <listitem>
369                     <para>
370                         This is for stiff problems. <literal>lsode</literal> solver of
371                         package ODEPACK is called and it uses the BDF method.
372                     </para>
373                 </listitem>
374             </varlistentry>
375             <varlistentry>
376                 <term>"rk":</term>
377                 <listitem>
378                     <para>Adaptive Runge-Kutta of order 4 (RK4) method.</para>
379                 </listitem>
380             </varlistentry>
381             <varlistentry>
382                 <term>"rkf":</term>
383                 <listitem>
384                     <para>
385                         The Shampine and Watts program based on Fehlberg's Runge-Kutta
386                         pair of order 4 and 5 (RKF45) method is used. This is for non-stiff
387                         and mildly stiff problems when derivative evaluations are
388                         inexpensive. This method should generally not be used when the user
389                         is demanding high accuracy.
390                     </para>
391                 </listitem>
392             </varlistentry>
393             <varlistentry>
394                 <term>"fix":</term>
395                 <listitem>
396                     <para>
397                         Same solver as <literal>"rkf"</literal>, but the user interface
398                         is very simple,
399                         i.e. only <varname>rtol</varname> and <varname>atol</varname>
400                         parameters can be passed to the solver. This is the simplest method
401                         to try.
402                     </para>
403                 </listitem>
404             </varlistentry>
405             <varlistentry>
406                 <term>"root":</term>
407                 <listitem>
408                     <para>
409                         ODE solver with rootfinding capabilities. The
410                         <literal>lsodar</literal> solver of package ODEPACK is used. It is a
411                         variant of the <literal>lsoda</literal> solver where it finds the
412                         roots of a given vector function. See help on <link linkend="ode_root">ode_root</link> for more
413                         details.
414                     </para>
415                 </listitem>
416             </varlistentry>
417             <varlistentry>
418                 <term>"discrete":</term>
419                 <listitem>
420                     <para>
421                         Discrete time simulation. See help on
422                         <link linkend="ode_discrete">ode_discrete</link> for more
423                         details.
424                     </para>
425                 </listitem>
426             </varlistentry>
427         </variablelist>
428     </refsection>
429     <refsection>
430         <title>Examples</title>
431         <para>
432             In the following example, we solve the Ordinary Differential Equation
433             <literal>dy/dt=y^2-y*sin(t)+cos(t)</literal> with the initial
434             condition <literal>y(0)=0</literal>.
435             We use the default solver.
436         </para>
437         <programlisting role="example"><![CDATA[
438 function ydot=f(t,y)
439         ydot=y^2-y*sin(t)+cos(t)
440 endfunction
441 y0=0;
442 t0=0;
443 t=0:0.1:%pi;
444 y = ode(y0,t0,t,f);
445 plot(t,y)
446  ]]></programlisting>
447         <scilab:image>
448             function ydot=f(t,y)
449             ydot=y^2-y*sin(t)+cos(t)
450             endfunction
451             y0=0;
452             t0=0;
453             t=0:0.1:%pi;
454             y = ode(y0,t0,t,f);
455             plot(t,y)
456         </scilab:image>
457         <para>
458             In the following example, we solve the equation <literal>dy/dt=A*y</literal>.
459             The exact solution is <literal>y(t)=expm(A*t)*y(0)</literal>, where
460             <literal>expm</literal> is the matrix exponential.
461             The unknown is the 2-by-1 matrix y(t).
462         </para>
463         <programlisting role="example"><![CDATA[
464 function ydot=f(t,y)
465     ydot=A*y
466 endfunction
467 function J=Jacobian(t,y)
468     J=A
469 endfunction
470 A=[10,0;0,-1];
471 y0=[0;1];
472 t0=0;
473 t=1;
474 ode("stiff",y0,t0,t,f,Jacobian)
475 // Compare with exact solution:
476 expm(A*t)*y0
477  ]]></programlisting>
478         <para>
479             In the following example, we solve the ODE <literal>dx/dt = A x(t) + B u(t)</literal>
480             with <literal>u(t)=sin(omega*t)</literal>.
481             Notice the extra arguments of <varname>f</varname>:
482             <literal>A</literal>, <literal>u</literal>, <literal>B</literal>,
483             <literal>omega</literal> are passed to function <varname>f</varname> in a list.
484         </para>
485         <programlisting role="example"><![CDATA[
486 function xdot=linear(t,x,A,u,B,omega)
487     xdot=A*x+B*u(t,omega)
488 endfunction
489 function ut=u(t,omega)
490     ut=sin(omega*t)
491 endfunction
492 A=[1 1;0 2];
493 B=[1;1];
494 omega=5;
495 y0=[1;0];
496 t0=0;
497 t=[0.1,0.2,0.5,1];
498 ode(y0,t0,t,list(linear,A,u,B,omega))
499  ]]></programlisting>
500         <para>
501             In the following example, we solve the Riccati differential equation
502             <literal>dX/dt=A'*X + X*A - X'*B*X + C</literal> where initial
503             condition <literal>X(0)</literal> is the 2-by-2 identity matrix.
504         </para>
505         <programlisting role="example"><![CDATA[
506 function Xdot=ric(t,X,A,B,C)
507     Xdot=A'*X+X*A-X'*B*X+C
508 endfunction
509 A=[1,1;0,2];
510 B=[1,0;0,1];
511 C=[1,0;0,1];
512 y0=eye(A);
513 t0=0;
514 t=0:0.1:%pi;
515 X = ode(y0,t0,t,list(ric,A,B,C))
516  ]]></programlisting>
517         <para>
518             In the following example, we solve the differential equation
519             <literal>dY/dt=A*Y</literal> where the unknown <literal>Y(t)</literal>
520             is a 2-by-2 matrix.
521             The exact solution is <literal>Y(t)=expm(A*t)</literal>, where
522             <literal>expm</literal> is the matrix exponential.
523         </para>
524         <programlisting role="example"><![CDATA[
525 function ydot=f(t,y,A)
526     ydot=A*y;
527 endfunction
528 A=[1,1;0,2];
529 y0=eye(A);
530 t0=0;
531 t=1;
532 ode(y0,t0,t,list(f,A))
533 // Compare with the exact solution:
534 expm(A*t)
535 ode("adams",y0,t0,t,f)
536  ]]></programlisting>
537     </refsection>
538     <refsection>
539         <title>With a compiler</title>
540         <para>
541             The following example requires a C compiler.
542         </para>
543         <programlisting role="example"><![CDATA[
544 // ---------- Simple one dimension ODE (C coded external)
545 ccode=['#include <math.h>'
546        'void myode(int *n,double *t,double *y,double *ydot)'
547        '{'
548        '  ydot[0]=y[0]*y[0]-y[0]*sin(*t)+cos(*t);'
549        '}']
550 mputl(ccode,TMPDIR+'/myode.c') //create the C file
551 // Compile
552 ilib_for_link('myode','myode.c',[],'c',[],TMPDIR+'/loader.sce');
553 exec(TMPDIR+'/loader.sce') //incremental linking
554 y0=0;
555 t0=0;
556 t=0:0.1:%pi;
557 y = ode(y0,t0,t,'myode');
558
559  ]]></programlisting>
560     </refsection>
561     <refsection role="see also">
562         <title>See Also</title>
563         <simplelist type="inline">
564             <member>
565                 <link linkend="ode_discrete">ode_discrete</link>
566             </member>
567             <member>
568                 <link linkend="ode_root">ode_root</link>
569             </member>
570             <member>
571                 <link linkend="dassl">dassl</link>
572             </member>
573             <member>
574                 <link linkend="impl">impl</link>
575             </member>
576             <member>
577                 <link linkend="odedc">odedc</link>
578             </member>
579             <member>
580                 <link linkend="odeoptions">odeoptions</link>
581             </member>
582             <member>
583                 <link linkend="csim">csim</link>
584             </member>
585             <member>
586                 <link linkend="ltitr">ltitr</link>
587             </member>
588             <member>
589                 <link linkend="rtitr">rtitr</link>
590             </member>
591             <member>
592                 <link linkend="intg">intg</link>
593             </member>
594         </simplelist>
595     </refsection>
596     <refsection>
597         <title>Bibliography</title>
598         <para>
599             Alan C. Hindmarsh, "lsode and lsodi, two new initial value ordinary
600             differential equation solvers", ACM-Signum newsletter, vol. 15, no. 4
601             (1980), pp. 10-11.
602         </para>
603     </refsection>
604     <refsection>
605         <title>Used Functions</title>
606         <para>
607             The associated routines can be found in SCI/modules/differential_equations/src/fortran directory:
608             lsode.f lsoda.f lsodar.f
609         </para>
610     </refsection>
611 </refentry>