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                         <literal>y</literal>, 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                         <literal>y</literal>, 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 <literal>f</literal>.
  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             <literal>ode</literal> 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 <literal>ode</literal> for
 166             standard explicit ODE systems.
 167         </para>
 168         <para>
 169             The simplest call of <literal>ode</literal> 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 <literal>y</literal> is computed and <literal>y</literal> is
 174             matrix of solution vectors
 175             <literal>y=[y(t(1)),y(t(2)),...]</literal>.
 176         </para>
 177         <para>
 178             The input argument <literal>f</literal> 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 <literal>f</literal> 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                     <literal>y</literal> 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 <literal>f</literal> 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, <literal>y</literal> 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, <literal>y</literal> 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 <literal>f</literal> needs extra
 237                     arguments.
 238                     In this case, we can use the following feature.
 239                     The <literal>f</literal> 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 <literal>f</literal> 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 <literal>ode</literal> 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 <literal>rtol</literal> and <literal>atol</literal> 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 <literal>rtol</literal> and/or
 266             <literal>atol</literal> is a constant <literal>rtol(i)</literal>
 267             and/or <literal>atol(i)</literal> are set to this constant value.
 268             Default values for <literal>rtol</literal> and <literal>atol</literal>
 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 <literal>jac</literal>.
 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 <literal>jac</literal> 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 <literal>y</literal> 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 <literal>f</literal>.
 290                 </para>
 291             </listitem>
 292             <listitem>
 293                 <para>
 294                     If <literal>jac</literal> 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 <literal>jac</literal> is a list the same conventions as for
 319                     <literal>f</literal> apply.
 320                 </para>
 321             </listitem>
 322         </itemizedlist>
 323         <para>
 324             Optional arguments <literal>w</literal> and
 325             <literal>iw</literal> 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 <literal>ode</literal> 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             <literal>type</literal> 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 <literal>rtol</literal> and <literal>atol</literal>
 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         <para>
 448             In the following example, we solve the equation <literal>dy/dt=A*y</literal>.
 449             The exact solution is <literal>y(t)=expm(A*t)*y(0)</literal>, where
 450             <literal>expm</literal> is the matrix exponential.
 451             The unknown is the 2-by-1 matrix y(t).
 452         </para>
 453         <programlisting role="example"><![CDATA[
 454 function ydot=f(t,y)
 455     ydot=A*y
 456 endfunction
 457 function J=Jacobian(t,y)
 458     J=A
 459 endfunction
 460 A=[10,0;0,-1];
 461 y0=[0;1];
 462 t0=0;
 463 t=1;
 464 ode("stiff",y0,t0,t,f,Jacobian)
 465 // Compare with exact solution:
 466 expm(A*t)*y0
 467  ]]></programlisting>
 468         <para>
 469             In the following example, we solve the ODE <literal>dx/dt = A x(t) + B u(t)</literal>
 470             with <literal>u(t)=sin(omega*t)</literal>.
 471             Notice the extra arguments of <literal>f</literal>:
 472             <literal>A</literal>, <literal>u</literal>, <literal>B</literal>,
 473             <literal>omega</literal> are passed to function <literal>f</literal> in a list.
 474         </para>
 475         <programlisting role="example"><![CDATA[
 476 function xdot=linear(t,x,A,u,B,omega)
 477     xdot=A*x+B*u(t,omega)
 478 endfunction
 479 function ut=u(t,omega)
 480     ut=sin(omega*t)
 481 endfunction
 482 A=[1 1;0 2];
 483 B=[1;1];
 484 omega=5;
 485 y0=[1;0];
 486 t0=0;
 487 t=[0.1,0.2,0.5,1];
 488 ode(y0,t0,t,list(linear,A,u,B,omega))
 489  ]]></programlisting>
 490         <para>
 491             In the following example, we solve the Riccati differential equation
 492             <literal>dX/dt=A'*X + X*A - X'*B*X + C</literal> where initial
 493             condition <literal>X(0)</literal> is the 2-by-2 identity matrix.
 494         </para>
 495         <programlisting role="example"><![CDATA[
 496 function Xdot=ric(t,X,A,B,C)
 497     Xdot=A'*X+X*A-X'*B*X+C
 498 endfunction
 499 A=[1,1;0,2];
 500 B=[1,0;0,1];
 501 C=[1,0;0,1];
 502 y0=eye(A);
 503 t0=0;
 504 t=0:0.1:%pi;
 505 X=ode(y0,t0,t,list(ric,A,B,C))
 506  ]]></programlisting>
 507         <para>
 508             In the following example, we solve the differential equation
 509             <literal>dY/dt=A*Y</literal> where the unknown <literal>Y(t)</literal>
 510             is a 2-by-2 matrix.
 511             The exact solution is <literal>Y(t)=expm(A*t)</literal>, where
 512             <literal>expm</literal> is the matrix exponential.
 513         </para>
 514         <programlisting role="example"><![CDATA[
 515 function ydot=f(t,y,A)
 516     ydot=A*y;
 517 endfunction
 518 A=[1,1;0,2];
 519 y0=eye(A);
 520 t0=0;
 521 t=1;
 522 ode(y0,t0,t,list(f,A))
 523 // Compare with the exact solution:
 524 expm(A*t)
 525 ode("adams",y0,t0,t,f)
 526  ]]></programlisting>
 527     </refsection>
 528     <refsection>
 529         <title>With a compiler</title>
 530         <para>
 531             The following example requires a C compiler.
 532         </para>
 533         <programlisting role="example"><![CDATA[
 534 // ---------- Simple one dimension ODE (C coded external)
 535 ccode=['#include <math.h>'
 536        'void myode(int *n,double *t,double *y,double *ydot)'
 537        '{'
 538        '  ydot[0]=y[0]*y[0]-y[0]*sin(*t)+cos(*t);'
 539        '}']
 540 mputl(ccode,TMPDIR+'/myode.c') //create the C file
 541 // Compile
 542 ilib_for_link('myode','myode.c',[],'c',TMPDIR+'/Makefile',TMPDIR+'/loader.sce');
 543 exec(TMPDIR+'/loader.sce') //incremental linking
 544 y0=0;
 545 t0=0;
 546 t=0:0.1:%pi;
 547 y=ode(y0,t0,t,'myode');
 548
 549  ]]></programlisting>
 550     </refsection>
 551     <refsection role="see also">
 552         <title>See Also</title>
 553         <simplelist type="inline">
 554             <member>
 555                 <link linkend="ode_discrete">ode_discrete</link>
 556             </member>
 557             <member>
 558                 <link linkend="ode_root">ode_root</link>
 559             </member>
 560             <member>
 561                 <link linkend="dassl">dassl</link>
 562             </member>
 563             <member>
 564                 <link linkend="impl">impl</link>
 565             </member>
 566             <member>
 567                 <link linkend="odedc">odedc</link>
 568             </member>
 569             <member>
 570                 <link linkend="odeoptions">odeoptions</link>
 571             </member>
 572             <member>
 573                 <link linkend="csim">csim</link>
 574             </member>
 575             <member>
 576                 <link linkend="ltitr">ltitr</link>
 577             </member>
 578             <member>
 579                 <link linkend="rtitr">rtitr</link>
 580             </member>
 581             <member>
 582                 <link linkend="intg">intg</link>
 583             </member>
 584         </simplelist>
 585     </refsection>
 586     <refsection>
 587         <title>Bibliography</title>
 588         <para>
 589             Alan C. Hindmarsh, "lsode and lsodi, two new initial value ordinary
 590             differential equation solvers", ACM-Signum newsletter, vol. 15, no. 4
 591             (1980), pp. 10-11.
 592         </para>
 593     </refsection>
 594     <refsection>
 595         <title>Used Functions</title>
 596         <para>
 597             The associated routines can be found in SCI/modules/differential_equations/src/fortran directory:
 598             lsode.f lsoda.f lsodar.f
 599         </para>
 600     </refsection>
 601 </refentry>