ode - ordinary differential equation solver
In this help we only describe the use of ode for standard explicit ODE systems.
The simplest call of ode is: y=ode(y0,t0,t,f) where y0 is the vector of initial conditions, t0 is the initial time, t is the vector of times at which the solution y is computed and y is matrix of solution vectors y=[y(t(1)),y(t(2)),...] .
The input argument f defines the RHS of the first order differential equation: dy/dt=f(t,y). It is an external i.e. a function with specified syntax, or the name of a Fortran subroutine or a C function (character string) with specified calling sequence or a list:
The Fortran routine must have the following calling sequence: fex(n,t,y,ydot) , with n an integer, t a double precision scalar, y and ydot double precision vectors.
The C function must have the following prototype: fex(int *n,double *t,double *y,double *ydot)
t is the time, y the state and ydot the state derivative (dy/dt)This external can be build in a OS independant way using ilib_for_link and dynamically linked to Scilab by the link function.
This syntax allows to use parameters as the arguments of realf .
The function f can return a p x q matrix instead of a vector. With this matrix notation, we solve the n=p+q ODE's system dY/dt=F(t,Y) where Y is a p x q matrix. Then initial conditions, Y0 , must also be a p x q matrix and the result of ode is the p x q(T+1) matrix [Y(t_0),Y(t_1),...,Y(t_T)] .
and integration is carried out as far as this error is small for all components of the state. If rtol and/or atol is a constant rtol(i) and/or atol(i) are set to this constant value. Default values for rtol and atol are respectively rtol=1.d-5 and atol=1.d-7 for most solvers and rtol=1.d-3 and atol=1.d-4 for "rfk" and "fix" .
For stiff problems, it is better to give the Jacobian of the RHS function as the optional argument jac . It is an external i.e. a function with specified syntax, or the name of a Fortran subroutine or a C function (character string) with specified calling sequence or a list.
If jac is a function the syntax should be J=jac(t,y)
where t is a real scalar (time) and y a real vector (state). The result matrix J must evaluate to df/dx i.e. J(k,i) = dfk/dxi with fk = kth component of f.
If jac is a character string it refers to the name of a Fortran subroutine or a C function, with the following calling sequence:
Fortran case:subroutine fex(n,t,y,ml,mu,J,nrpd) integer n,ml,mu,nrpd double precision t,y(*),J(*)C case:
void fex(int *n,double *t,double *y,int *ml,int *mu,double *J,int *nrpd,)
jac(n,t,y,ml,mu,J,nrpd) . In most cases you have not to refer ml , mu and nrpd .
If jac is a list the same conventions as for f apply.
More options can be given to ODEPACK solvers by using %ODEOPTIONS variable. See odeoptions .
// ---------- Simple one dimension ODE (Scilab function external)
// dy/dt=y^2-y sin(t)+cos(t), y(0)=0
function ydot=f(t,y),ydot=y^2-y*sin(t)+cos(t),endfunction
y0=0;t0=0;t=0:0.1:%pi;
y=ode(y0,t0,t,f)
plot(t,y)
// ---------- Simple one dimension ODE (C coded external)
ccode=['#include <math.h>'
'void myode(int *n,double *t,double *y,double *ydot)'
'{'
' ydot[0]=y[0]*y[0]-y[0]*sin(*t)+cos(*t);'
'}']
mputl(ccode,TMPDIR+'/myode.c') //create the C file
ilib_for_link('myode','myode.o',[],'c',TMPDIR+'/Makefile',TMPDIR+'/loader.sce');//compile
exec(TMPDIR+'/loader.sce') //incremental linking
y0=0;t0=0;t=0:0.1:%pi;
y=ode(y0,t0,t,'myode');
// ---------- Simulation of dx/dt = A x(t) + B u(t) with u(t)=sin(omega*t),
// x0=[1;0]
// solution x(t) desired at t=0.1, 0.2, 0.5 ,1.
// A and u function are passed to RHS function in a list.
// B and omega are passed as global variables
function xdot=linear(t,x,A,u),xdot=A*x+B*u(t),endfunction
function ut=u(t),ut=sin(omega*t),endfunction
A=[1 1;0 2];B=[1;1];omega=5;
ode([1;0],0,[0.1,0.2,0.5,1],list(linear,A,u))
// ---------- Matrix notation Integration of the Riccati differential equation
// Xdot=A'*X + X*A - X'*B*X + C , X(0)=Identity
// Solution at t=[1,2]
function Xdot=ric(t,X),Xdot=A'*X+X*A-X'*B*X+C,endfunction
A=[1,1;0,2]; B=[1,0;0,1]; C=[1,0;0,1];
t0=0;t=0:0.1:%pi;
X=ode(eye(A),0,t,ric)
//
// ---------- Matrix notation, Computation of exp(A)
A=[1,1;0,2];
function xdot=f(t,x),xdot=A*x;,endfunction
ode(eye(A),0,1,f)
ode("adams",eye(A),0,1,f)
// ---------- Matrix notation, Computation of exp(A) with stiff matrix, Jacobian given
A=[10,0;0,-1];
function xdot=f(t,x),xdot=A*x,endfunction
function J=Jacobian(t,y),J=A,endfunction
ode("stiff",[0;1],0,1,f,Jacobian)
ode_discrete , ode_root , dassl , impl , odedc , odeoptions , csim , ltitr , rtitr ,
Alan C. Hindmarsh, lsode and lsodi, two new initial value ordinary differential equation solvers, acm-signum newsletter, vol. 15, no. 4 (1980), pp. 10-11.
The associated routines can be found in routines/integ directory :
lsode.f lsoda.f lsodar.f