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>