[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