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

[Karline Soetaert]

Hi Maciek,
I think it is because you request such detailed output. Try:

times <- seq(0, 300, by = 1)
rather than 
times <- seq(0, 300, by = 0.01)


Karline Soetaert
-----Original Message-----
From: R-sig-dynamic-models [mailto:r-sig-dynamic-models-bounces at r-project.org] On Behalf Of Maciek Jacek Swat
Sent: woensdag 22 november 2017 10:40
To: r-sig-dynamic-models at r-project.org
Subject: [R-sig-dyn-mod] lsoda (deSolve) order of magnitude slower than ode45

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