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