The different ways of coding user functions used by IDA
To solve a problem with IDA different functions have to be provided. For example the residual function of the dae has to be computed by the user and is the only mandatory function. Providing the Jacobian always help the solver, for example when the ode is stiff. Two other optional functions allow to handle events and intermediate callbacks (at each user prescribed or internal solver step). The usual way is to write functions in the Scilab language but if speed is a concern, entrypoints of dynamically linked shared libraries (DLL, built from C,C++ of Fortran) should be used instead, as at least one order of magnitude in speed is typically achieved that way (examples below use C). In that case all arguments are given via their address (Fortran convention).
Extra parameters can be passed to user Scilab functions or entrypoints by replacing the corresponding argument or option by a list, where the first element is the usual argument (a function identifier or a string with the entrypoint name) and the subsequent elements are the parameters to be passed after the mandatory arguments. Note that DLL entrypoints accept only one extra parameter as a double array.
In the call ida(f,tspan,y0,yp0)
the first argument f
can be:
a function with prototype eq=f(t,y,yp,...)
function eq=f(t, y, yp) eq = yp+y-1 end function eq=fp(t, y, yp, a, b) eq = yp-(a*y+b); end [t,y,yp] = ida(f,[0 10],0,1); [t,y,yp] = ida(list(fp,-1,1),[0 10],0,1); | ![]() | ![]() |
a string, for example "dynres"
defining the name of a DLL entrypoint
with C prototype
void dynrhs(int *n,double *t,double y[],double dydt[],double res[])
The following example shows the complete workflow up to the final solver call:
rescode=["void dynres(int *n,double *t,double y[],double yp[],double res[])" "{" "res[0] = yp[0]+y[0]-1;" "}" "void dynrespar(int *n,double *t,double y[],double yp[],double res[],double par[])" "{" "res[0] = yp[0]-(par[0]*y[0]+par[1]);" "}"] cd(TMPDIR) mputl(rescode,"dynres.c") //create the C file ilib_for_link(["dynres","dynrespar"],"dynres.c",[],"c");//compile exec("loader.sce",-1) [t,y,yp] = ida("dynres",[0 10],0,1); [t,y,yp] = ida(list("dynrespar",[-1,1]),[0 10],0,1); | ![]() | ![]() |
In IDA Newton iterations the Jacobian of the residual system is by default
approximated by finite differences. Iterations can be accelerated and precision can be improved by giving the true
Jacobian of the residual with respect to y and yp with the jacobian
option. In the call
ida(f,tspan,y0,yp0,jacobian=jacyyp)
, jacyyp
can be:
a function with prototype [jacy,jacyp]=jacyyp(t,y,yp,...)
where jacy,jacyp
are
the Jacobian of the residual function with respect to y and yp, respectively. For example for the SIR system:
function res=sir(t, y, yp, bet, gam) res=[yp(1)+bet*y(1)*y(2) yp(2)-bet*y(1)*y(2)+gam*y(2) y(1)+y(2)+y(3)-1]; end function [jacy, jacyp]=sir_jac(t, y, yp, bet, gam) jacy = [ bet*y(2) bet*y(1) 0 -bet*y(2) -bet*y(1)+gam 0 1 1 1]; jacyp = diag([1 1 0]); end gam=1/40; bet=0.2; y0=[0.999;0.001;0]; yp0 = [-bet*y0(1)*y0(2);+bet*y0(1)*y0(2)-gam*y0(2);gam*y0(2)]; sol = ida(list(sir,bet,gam),[0 200],y0,yp0,jacobian=list(sir_jac,bet,gam)); | ![]() | ![]() |
a cell containing the constant Jacobians with respect to y and yp, which allows to easily consider the case of a linear DAE, e.g.
a string giving the name of a DLL entrypoint with C prototype
void dynjac(int *n,double *t,double *cj,double y[],double ydot[],double res[],double J[])
The current residual f(t,y,yp)
is passed by convenience as it may be usefull in the
computation. The output J
is the sum of the Jacobian w.r.t. y and *c
times the Jacobian w.r.t. yp, stored column-wise. Here is a full example (Robertson DAE):
code=['void chemres(int *n, double *t, double y[], double yd[], double r[])' ' {' ' r[0] = yd[0]+0.04*y[0]-1.0e4*y[1]*y[2];' ' r[1] = yd[1]-0.04*y[0]+1.0e4*y[1]*y[2]+3.0e7*y[1]*y[1];' ' r[2] = y[0]+y[1]+y[2]-1;' ' }' 'void chemjac(int *n, double *t, double *cj, double y[], double yd[], double r[], double j[])' ' {' ' /* first column*/' ' j[0] = 0.04+*cj;' ' j[1] = -0.04;' ' j[2] = 1.0;' ' /* second column*/' ' j[3] = -1.0e4*y[2];' ' j[4] = +1.0e4*y[2]+2*3.0e7*y[1]+*cj;' ' j[5] = 1.0;' ' /* third column*/' ' j[6] = -1.0e4*y[1];' ' j[7] = +1.0e4*y[1];' ' j[8] = 1.0;' ' }']; y0 = [1-1e-3; 0; 1e-3]; yp0 = [-0.0400; 0.0400; 0]; tspan = [0 4e6]; mputl(code, 'code.c'); ilib_for_link(['chemres','chemjac'],'code.c',[],'c'); exec('loader.sce',-1); [t,y,yp] = ida("chemres",tspan,y0,yp0,jacobian="chemjac"); y(2,:)=y(2,:)*20000; semilogx(t(2:$),y(:,2:$)) | ![]() | ![]() |
In the call ida(f,tspan,y0,yp0,events=g)
, g
can be:
a function with prototype [eq,term,dir] = g(t,y,yp,...)
, as described in the Events section of the IDA main page.
a string giving the name of a DLL entrypoint with C prototype
void dynevent(int *n,double *t,double y[], double yp[], int *ng, double g[], int term[], int dir[])
When initializing the solver the entrypoint is called with *ng == 0
and *ng
must be set to the number of events before exiting. When g==NULL
the two arrays term
and dir
must be filled to define the termination
conditions and events directions (see the relevant values on the Events section of the
CVODE main page). Otherwise the entrypoint must return the event vector in the array g
.
See the example below:
CVODE can call a Scilab function or DLL entryoint after every successfull internal or user prescribed step. In the
call ida(f,tspan,y0,yp0,intcb=callback)
, callback
can be:
a function with prototype term = callback(t,y,yp,flag,stats,...)
where flag
can take the values
"init"
, "step"
or "done"
and stats
is a structure gathering the solver statistics. The function has to return %f
to continue
integration whereas returning %t
will stop the solver.
a string giving the name of a DLL entrypoint with prototype
int dyncb(int *n,double *t,double y[],double yp[],int* flag)
The entrypoint has to return 0 or 1, for normal continuation or solver termination, respectively. The flag argument takes values -1,0 or 1 for "init", "step" and "done" solver phases.