[R-sig-dyn-mod] lsodes with jacobian

Daniel Kaschek daniel.kaschek at physik.uni-freiburg.de
Sat Jun 28 16:31:29 CEST 2014


Hi again,

by looking at call_losda.c in the source code of the deSolve package, I
found the necessary call. You can find my fist example at the end of the
mail.

Next, I tried this out for an ODE with 600 equations. Unfortunately, it
takes 60 times longer with jacobian than without :-(
I suppose, this has to do with the sparsity specifications. Or is it the
"if"-lines in the jacobian? 

* What would be the correct specification of sparsetype, nnz and inz in
case of an explicitly given jacvec function?
* Knowing the zeros of my jacobian, does it make sens to write down the
jacvec function explicitly or is it sufficient to specify nnz and inz?

Cheers,
Daniel

#include <R.h> 
 #include <math.h> 

static double parms[6];
static double forc[2];

#define rGrow parms[0] 
 #define K parms[1] 
 #define rIng parms[2] 
 #define rMort parms[3] 
 #define y0_0 parms[4] 
 #define y1_0 parms[5] 
 #define assEff forc[0] 
 #define assEffDot forc[1] 

void initmod(void (* odeparms)(int *, double *)) {
	 int N=6;
	 odeparms(&N, parms);
}

void initforc(void (* odeforcs)(int *, double *)) {
	 int N=2;
	 odeforcs(&N, forc);
}

/** Derivatives (ODE system) **/
void derivs (int *n, double *t, double *y, double *ydot, double *RPAR,
int *IPAR) {


	 ydot[0] = rGrow*y[0]*(1-y[0]/K)-rIng*y[0]*y[1];
 	 ydot[1] = rIng*y[0]*y[1]*assEff-rMort*y[1];

}

/** Jacobian vector of the ODE system **/
void jacvec (int *neq, double *t, double *y, int *j, int *ian, int *jan,
double *pdj, double *yout, int *iout) {

	 int i;
for(i=0; i<*neq; i++) pdj[i] = 0.;
	 if(*j == 0 ) {
	 pdj[0] = rGrow*(1-y[0]/K)-rGrow*y[0]*(1/K)-rIng*y[1];
 	 pdj[1] = rIng*y[1]*assEff;
	 }
	 else if(*j == 1 ) {
	 pdj[0] = -(rIng*y[0]);
 	 pdj[1] = rIng*y[0]*assEff-rMort;
	 }

}



More information about the R-sig-dynamic-models mailing list