1 <?xml version="1.0" encoding="UTF-8"?>
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
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
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">
17 <refname>ode</refname>
18 <refpurpose>ordinary differential equation solver</refpurpose>
21 <title>Calling Sequence</title>
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)
30 <title>Arguments</title>
35 <para>a real vector or matrix, the initial conditions.</para>
42 a real scalar, the initial time.
50 a real vector, the times at which the solution is computed.
58 a function, external, string or list, the right hand side of
59 the differential equation.
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>.
79 a real constant or real vector of the same size as
80 <literal>y</literal>, the relative tolerance.
88 a real constant or real vectors of the same size as
89 <literal>y</literal>, the absolute tolerance.
97 a function, external, string or list, the Jacobian of the
98 function <literal>f</literal>.
105 <para>real vectors. (INPUT/OUTPUT)</para>
111 <para>an integer.</para>
117 <para>an external (function or character string or list).</para>
123 <para>an integer (initial time).</para>
129 <para>an integer vector.</para>
135 <para>a real vector or matrix. (OUTPUT)</para>
141 <para>a real vector. (OUTPUT)</para>
147 <title>Description</title>
149 <literal>ode</literal> solves explicit Ordinary Different Equations defined by:
155 \dfrac{dy}{dt} = f(t,y)\\
162 It is an interface to various solvers, in particular to ODEPACK.
165 In this help, we only describe the use of <literal>ode</literal> for
166 standard explicit ODE systems.
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>.
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.
184 If <literal>f</literal> is a Scilab function, its calling sequence
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
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.
205 The Fortran routine must have the header:
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
217 The C function must have the header:
220 fex(int *n,double *t,double *y,double *ydot)
223 where <literal>t</literal> is the time, <literal>y</literal> the
224 state and <literal>ydot</literal> is the state derivative
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>
236 It may happen that the simulator <literal>f</literal> needs extra
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>.
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>.
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>.
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.
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>.
294 If <literal>jac</literal> is a character string it refers to the
295 name of a Fortran subroutine or a C function.
298 The Fortran routine must have the header:
301 subroutine fex(n,t,y,ml,mu,J,nrpd)
303 double precision t,y(*),J(*)
306 The C function must have the header:
309 void fex(int *n,double *t,double *y,int *ml,int *mu,double *J,int *nrpd,)
312 In most cases you have not to refer <literal>ml</literal>, <literal>mu</literal> and
313 <literal>nrpd</literal>.
318 If <literal>jac</literal> is a list the same conventions as for
319 <literal>f</literal> apply.
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
332 More options can be given to ODEPACK solvers by using
333 <literal>%ODEOPTIONS</literal> variable. See <link linkend="odeoptions">odeoptions</link>.
337 <title>The solvers</title>
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:
345 <term><not given>:</term>
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
358 <term>"adams":</term>
361 This is for nonstiff problems. <literal>lsode</literal> solver
362 of package ODEPACK is called and it uses the Adams method.
367 <term>"stiff":</term>
370 This is for stiff problems. <literal>lsode</literal> solver of
371 package ODEPACK is called and it uses the BDF method.
378 <para>Adaptive Runge-Kutta of order 4 (RK4) method.</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.
397 Same solver as <literal>"rkf"</literal>, but the user interface
399 i.e. only <literal>rtol</literal> and <literal>atol</literal>
400 parameters can be passed to the solver. This is the simplest method
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
418 <term>"discrete":</term>
421 Discrete time simulation. See help on
422 <link linkend="ode_discrete">ode_discrete</link> for more
430 <title>Examples</title>
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.
437 <programlisting role="example"><![CDATA[
439 ydot=y^2-y*sin(t)+cos(t)
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).
453 <programlisting role="example"><![CDATA[
457 function J=Jacobian(t,y)
464 ode("stiff",y0,t0,t,f,Jacobian)
465 // Compare with exact solution:
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.
475 <programlisting role="example"><![CDATA[
476 function xdot=linear(t,x,A,u,B,omega)
477 xdot=A*x+B*u(t,omega)
479 function ut=u(t,omega)
488 ode(y0,t0,t,list(linear,A,u,B,omega))
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.
495 <programlisting role="example"><![CDATA[
496 function Xdot=ric(t,X,A,B,C)
497 Xdot=A'*X+X*A-X'*B*X+C
505 X=ode(y0,t0,t,list(ric,A,B,C))
508 In the following example, we solve the differential equation
509 <literal>dY/dt=A*Y</literal> where the unknown <literal>Y(t)</literal>
511 The exact solution is <literal>Y(t)=expm(A*t)</literal>, where
512 <literal>expm</literal> is the matrix exponential.
514 <programlisting role="example"><![CDATA[
515 function ydot=f(t,y,A)
522 ode(y0,t0,t,list(f,A))
523 // Compare with the exact solution:
525 ode("adams",y0,t0,t,f)
529 <title>With a compiler</title>
531 The following example requires a C compiler.
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)'
538 ' ydot[0]=y[0]*y[0]-y[0]*sin(*t)+cos(*t);'
540 mputl(ccode,TMPDIR+'/myode.c') //create the C file
542 ilib_for_link('myode','myode.c',[],'c',TMPDIR+'/Makefile',TMPDIR+'/loader.sce');
543 exec(TMPDIR+'/loader.sce') //incremental linking
547 y=ode(y0,t0,t,'myode');
551 <refsection role="see also">
552 <title>See Also</title>
553 <simplelist type="inline">
555 <link linkend="ode_discrete">ode_discrete</link>
558 <link linkend="ode_root">ode_root</link>
561 <link linkend="dassl">dassl</link>
564 <link linkend="impl">impl</link>
567 <link linkend="odedc">odedc</link>
570 <link linkend="odeoptions">odeoptions</link>
573 <link linkend="csim">csim</link>
576 <link linkend="ltitr">ltitr</link>
579 <link linkend="rtitr">rtitr</link>
582 <link linkend="intg">intg</link>
587 <title>Bibliography</title>
589 Alan C. Hindmarsh, "lsode and lsodi, two new initial value ordinary
590 differential equation solvers", ACM-Signum newsletter, vol. 15, no. 4
595 <title>Used Functions</title>
597 The associated routines can be found in SCI/modules/differential_equations/src/fortran directory:
598 lsode.f lsoda.f lsodar.f