Skip to content
Prev 379 / 696 Next

[R-sig-dyn-mod] forcings {deSolve}

Hi Daniel, 
 
Many thanks for your detailed answer. Combining the matrices for "a" and "b" as a list was what I was missing.
 
Thanks for your help, 
 
All the best, 
Tom
Hi Tom,

even though you have several forcings, you have only one function
"initforc". Consider the following example of dx/dt=-k1*x+k2*a+k3*b
where "a" and "b" are forcings. The corresponding C code would be


/** --------- Auto-generated code ------------ **/

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

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

#define k1 parms[0] 
#define k2 parms[1] 
#define k3 parms[2] 
#define y0_0 parms[3] 
#define a forc[0] 
#define b forc[1] 

void initmod(void (* odeparms)(int *, double *)) {
	 int N=4;
	 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] = -k1*y[0]+k2*a+k3*b;

}


Then you call ode() with func = "derivs", initfunc = "initmod", initforc
= "initforc" and forcings = myforclist where myforclist may be

$a
	 t x
[1,] 0 1
[2,] 1 1
[3,] 2 1
[4,] 3 1

$b
	 t x
[1,] 0 1
[2,] 1 2
[3,] 2 3
[4,] 3 4

i.e. a list of length two. Each list entry is a matrix with two columns,
time and value. The list should be ordered according to the appearance
of forcings in the C file.

Sometimes the forcing is known as an algebraic function, e.g. an
exponential function. In this case you can define the forcing directly
in deriv():

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

	double a = exp(-k4 * *t);
	double b = exp(-k5 * *t);
	ydot[0] = -k1*y[0]+k2*a+k3*b;

}

In this case add k4 and k5 to the parameters and remove a and b from the
forcings.

Hope this helps.

Cheers,
Daniel
On Mi, 2015-03-04 at 15:36 +0000, Tom Sumner wrote: