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>