[R-sig-dyn-mod] lsoda (deSolve) order of magnitude slower than ode45

[Maciek Jacek Swat]

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]]