IDA user functions

Introduction

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).

Passing extra parameters

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.

Residual function

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);

Jacobian

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.

function out=res(t, y, yp, A, B)
    out = A*yp-B*y;
end
A = [2 1; 1 1];
B = [0 1;-1 -1];
y0 = [1;0];
yp0 = A\(B*y0);
sol = ida(list(res,A,B),[0 10],y0,yp0);
solj = ida(list(res,A,B),[0 10],y0,yp0,jacobian={-B,A});

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:$))

Events

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:

Callback

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.

Report an issue
<< IDA	sci_sundials	CVODE, ARKODE and IDA options >>