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