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

[Maciek Jacek Swat]

Thank you Karline! Now, this reduces the simulation to 0.05 seconds, i.e.
twice as fast as ode45.
Best Regards, Maciej


On Wed, Nov 22, 2017 at 9:51 AM, Karline Soetaert <Karline.Soetaert at nioz.nl>
wrote:

> 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
>
> _______________________________________________
> R-sig-dynamic-models mailing list
> R-sig-dynamic-models at r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-dynamic-models
>

	[[alternative HTML version deleted]]