[R-sig-dyn-mod] Individual Based Model type with deSolve
Romain Pete
rpete_mp1 at ifremer.fr
Wed Jan 31 10:54:47 CET 2018
Hi Folks,
I’m trying to implement an IBM type model in an existing biogeochemical model (NPZD-like). What I got is an embedded (NOT nested function) DEB oyster model (in my NPZD) on which I’d like to use a for loop to solve DEB equations for multiple seeding. Finally, I’d like to get the output of DEB with the output of the NPZD:
MyModel <- function(t, state, parameters) {
with(as.list(c(state, parameters)), {
# all biogeochemical processes defined here
for( s in seeding){
run DEB equations
c(seeding1, seeding2, etc..) # 4 state variables DEB hence a matrix result with nrow=4 state variables, ncol= « s » seeding
}
# rate of change
dN = remineralisation + oyster excretion - uptake
dP = growth - Zoo_grazing - Oyster_grazing - death
dZ = growth - Oyster_grazing - death - egestion
dD = death + egestion - mineralization
dE = stock ###
dV = volume ### These 4 state variables are multiplied by number of seedings
dEr = reproduction ###
dEgo = gamete ###
list(c(
dN, dP, dZ, dD, dO, c(seeding)
))
})
}
Trouble here is that the ode function uses initialization vector of length = number of state variables and when using multiple seedings, I, somehow, add state variables, hence, I’m struggling…
Any thoughts would be much appreciated
Cheers
Romain
> Le 22 nov. 2017 à 10:39, Maciek Jacek Swat <maciej.swat at gmail.com> a écrit :
>
> hi All,
>
> I ran a simple ODE model (see PS for R and Matlab code) and get order of
> magnitude longer running times with lsoda (as implemented in deSolve) than
> with ode45 in Matlab. (lsoda takes 2.97 seconds, ode45 takes 0.12 seconds).
> Is this so or can it be optimized?
>
> For example, the ODEs in R are coded using the notation dA <- ..., dB <-
> ... and not the vectorized one, is it the reason for the longer simulation
> times?
>
> Any comments would be very appreciated,
> Best, Maciej
>
>
> PS:
> R code:
> rm(list=ls())
>
> Nanda <- function(t, state, parameters) {
> with(as.list(c(state, parameters)), {
>
> dB_CLL <- 1/body*(((r*B_CLL)*body) - ((d_B*B_CLL)*body) + ((b_B)*body)
> - ((d_BN*B_CLL*NK)*body) - ((d_BT*B_CLL*Tc)*body))
> dNK <- 1/body*(-((d_NB*NK*B_CLL)*body) + ((b_N)*body) -
> ((d_N*NK)*body))
> dTc <- 1/body*(-((d_TB*Tc*B_CLL)*body) - ((d_T*Tc)*body) + ((b_T)*body)
> + (k*a_TH*B_CLL^L/(s+B_CLL^L)*TH*Tc))
> dTH <- 1/body*((a_TH*B_CLL^L/(s+B_CLL^L)*TH) - ((d_TH*TH)*body) +
> ((b_TH)*body))
>
> list(c(
> dB_CLL,dNK,dTc,dTH
> ))
> })
> }
>
> parameters <- c(
> a_TH = 0.0001,
> L = 2,
> s = 1000,
> d_B = 0.0076,
> r = 0.0176,
> b_B = 68,
> B_CLL_0 = 100,
> NK_0 = 1000,
> Tc_0 = 100,
> TH_0 = 100,
> d_BN = 0.0001,
> d_BT = 0.026,
> d_NB = 0.0001,
> b_N = 3.98,
> d_N = 0.0159,
> d_TB = 0.001,
> d_T = 0.005,
> b_T = 2.6,
> k = 0.6,
> d_TH = 0.01,
> b_TH = 10.4,
> body = 1
> )
>
> state <- c(
> B_CLL = 100,
> NK = 1000,
> Tc = 100,
> TH = 100
> )
>
> t1=proc.time()
> times <- seq(0, 300, by = 0.01)
> out <- ode(y = state, times = times, func = Nanda, method="lsoda", parms =
> parameters)
> cat((proc.time()-t1)[3],'\n')
>
>
>
>
> Matlab code:
> clear all; close all; clc;
>
> yy0(1) = 100; % B_CLL
> yy0(2) = 1000; % NK
> yy0(3) = 100; % Tc
> yy0(4) = 100; % TH
>
> tic
> [T,Y] = ode45(@(t,y)NandaFCT(t,y),[0 300],yy0);
> toc
>
> figure
> semilogy(T,Y); grid on;
>
>
> function dydt = NandaFCT(t,y)
>
> a_TH = 0.0001;
> L = 2;
> s = 1000;
> d_B = 0.0076;
> r = 0.0176;
> b_B = 68;
> B_CLL_0 = 100;
> NK_0 = 1000;
> Tc_0 = 100;
> TH_0 = 100;
> d_BN = 0.0001;
> d_BT = 0.026;
> d_NB = 0.0001;
> b_N = 3.98;
> d_N = 0.0159;
> d_TB = 0.001;
> d_T = 0.005;
> b_T = 2.6;
> k = 0.6;
> d_TH = 0.01;
> b_TH = 10.4;
> body = 1;
>
> dydt = zeros(4,1);
>
> dydt(1)= 1/body*(((r*y(1))*body) - ((d_B*y(1))*body) + ((b_B)*body) -
> ((d_BN*y(1)*y(2))*body) - ((d_BT*y(1)*y(3))*body)); % B_CLL
> dydt(2)= 1/body*(-((d_NB*y(2)*y(1))*body) + ((b_N)*body) -
> ((d_N*y(2))*body)); % NK
> dydt(3)= 1/body*(-((d_TB*y(3)*y(1))*body) - ((d_T*y(3))*body) +
> ((b_T)*body) + (k*a_TH*y(1)^L/(s+y(1)^L)*y(4)*y(3))); % Tc
> dydt(4)= 1/body*((a_TH*y(1)^L/(s+y(1)^L)*y(4)) - ((d_TH*y(4))*body) +
> ((b_TH)*body)); % TH
>
> end
>
> [[alternative HTML version deleted]]
>
> _______________________________________________
> R-sig-dynamic-models mailing list
> R-sig-dynamic-models at r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-dynamic-models
More information about the R-sig-dynamic-models
mailing list