Karline,


Yes, you're right, the model is defined as a DAE system of index 1. I tried the solver radau providing the mass matrix and it worked: I obtained the same results computed by the matlab ode15s. Unfortunately, there is still a significant computational difference:

ode15s: 35s
radau: 1333s

I'm going to try also with daspk...

Thanks for your help!


Ettore Mosca, PhD

CNR-Institute for Biomedical Technologies (http://www.itb.cnr.it)
Via Fratelli Cervi, 93, 20090, Segrate (Milan), Italy
Tel: +39 0226422600
Fax: +39 0226422660
skype: ettore.mosca


On 12/30/2011 12:00 PM, r-sig-dynamic-models-request@r-project.org wrote:
> ----------------------------------------------------------------------
>
> Message: 1
> Date: Thu, 29 Dec 2011 16:07:06 +0100
> From:"ettore.mosca@itb.cnr.it"  <ettore.mosca@itb.cnr.it>
> To:r-sig-dynamic-models@r-project.org
> Subject: [R-sig-dyn-mod] lsodes VS ode15s
> Message-ID:<4EFC821A.5000407@itb.cnr.it>
> Content-Type: text/plain
>
> Hello!
>
> First of all congratulations for the development of this very useful R
> package!
>
> I work with biochemical systems represented with differential equations
> and I have just successfully ported one model (approximately 60
> differential equations) from matlab to R: using lsodes and setting the
> appropriate rtol and atol values I obtain (almost the) same solutions of
> matlab ode15s and octave lsode.
>
> However, I'm not having the same luck with another model (same
> cardinality but more complex equations). While matlab ode15s solves it
> in 10s both octave lsode and R lsodes terminate the integration
> reporting various warnings (see below for R lsodes warnings) after 1 hour.
>
> I checked that initial conditions and parameters are the same; moreover,
> the function that solves the differential equations returns the "same"
> values of the matlab version: relative error = (y'_lsode - y'_matlab) /
> y'_matlab = 1e-6)
>
> Do you expect that lsodes cannot integrate some problems? I read, but I
> don't have expertise in this field, that while ode15s is based on the
> /numerical/  differentiation formulas, lsodes is based on the backward
> differentiation formula; maybe this is related with my issue?
>
> Or maybe I have to change some other integrator options?
>
> These are the matlab integrator options:
>
> odeset('RelTol',1e-3,'AbsTol',AbsTol_vec,'Vectorized','on',
> 'MassSingular','yes','Mass',M,'MStateDependence','none',
> 'MassConstant','yes','Events',@fevents);
>
> the two events implemented in the fevents function stop the integration
> when (i) the cpu time exceed a given threshold or (ii) a specific
> variable of the model is no longer biologically meaningful.
>
>
> thanks in advice for any suggestion,
>
>
> Ettore M.
>
>
> DLSODES- Warning..Internal T (=R1) and H (=R2) are
>         such that in the machine, T + H = T on the next step
>        (H = step size). Solver will continue anyway.
>         In above,  R1 =  0.2760657315701D-03   R2 =  0.1153501676498D-22
> DLSODES- Warning..Internal T (=R1) and H (=R2) are
>         such that in the machine, T + H = T on the next step
>        (H = step size). Solver will continue anyway.
>         In above,  R1 =  0.2760657315701D-03   R2 =  0.1153501676498D-22
> DLSODES- Warning..Internal T (=R1) and H (=R2) are
>         such that in the machine, T + H = T on the next step
>        (H = step size). Solver will continue anyway.
>         In above,  R1 =  0.2760657315701D-03   R2 =  0.3519053375292D-23
> DLSODES- Warning..Internal T (=R1) and H (=R2) are
>         such that in the machine, T + H = T on the next step
>        (H = step size). Solver will continue anyway.
>         In above,  R1 =  0.2760657315701D-03   R2 =  0.3519053375292D-23
> DLSODES- Warning..Internal T (=R1) and H (=R2) are
>         such that in the machine, T + H = T on the next step
>        (H = step size). Solver will continue anyway.
>         In above,  R1 =  0.2760657315701D-03   R2 =  0.8797633438231D-24
> DLSODES- Warning..Internal T (=R1) and H (=R2) are
>         such that in the machine, T + H = T on the next step
>        (H = step size). Solver will continue anyway.
>         In above,  R1 =  0.2760657315701D-03   R2 =  0.2199408359558D-24
> DLSODES- Warning..Internal T (=R1) and H (=R2) are
>         such that in the machine, T + H = T on the next step
>        (H = step size). Solver will continue anyway.
>         In above,  R1 =  0.2760657315701D-03   R2 =  0.2199408359558D-24
> DLSODES- Warning..Internal T (=R1) and H (=R2) are
>         such that in the machine, T + H = T on the next step
>        (H = step size). Solver will continue anyway.
>         In above,  R1 =  0.2760657315701D-03   R2 =  0.4158664893404D-24
> DLSODES- Warning..Internal T (=R1) and H (=R2) are
>         such that in the machine, T + H = T on the next step
>        (H = step size). Solver will continue anyway.
>         In above,  R1 =  0.2760657315701D-03   R2 =  0.1039666223351D-24
> DLSODES- Warning..Internal T (=R1) and H (=R2) are
>         such that in the machine, T + H = T on the next step
>        (H = step size). Solver will continue anyway.
>         In above,  R1 =  0.2760657315701D-03   R2 =  0.1039666223351D-24
> DLSODES- Above warning has been issued I1 times.
>        It will not be issued again for this problem.
>         In above message,  I1 =        10
>
>
> DLSODES- At current T (=R1), MXSTEP (=I1) steps
>         taken on this call before reaching TOUT
>         In above message,  I1 =      5000
>         In above message,  R1 =  0.2946239972790D-03
> Warning messages:
> 1: In lsodes(x0, DAEs, times = tspan, parms = p, atol = AbsTol_vec,  :
>     an excessive amount of work (>  maxsteps ) was done, but integration
> was not successful - increase maxsteps
> 2: In lsodes(x0, DAEs, times = tspan, parms = p, atol = AbsTol_vec,  :
>     Returning early. Results are accurate, as far as they go
>
>
>
> -- Ettore Mosca, PhD CNR-Institute for Biomedical Technologies 
> (http://www.itb.cnr.it) Via Fratelli Cervi, 93, 20090, Segrate 
> (Milan), Italy Tel: +39 0226422600 Fax: +39 0226422660 skype: 
> ettore.mosca [[alternative HTML version deleted]] 
> ------------------------------ Message: 2 Date: Thu, 29 Dec 2011 
> 17:02:42 +0100 From: "Soetaert, Karline" <K.Soetaert@nioo.knaw.nl> To: 
> "Special Interest Group for Dynamic Simulation Models in R" 
> <r-sig-dynamic-models@r-project.org> Subject: Re: [R-sig-dyn-mod] 
> lsodes VS ode15s Message-ID: 
> <65F6E1EC64DCA6489800C09A2007FC6E04090D28@cememail1.nioo.int> 
> Content-Type: text/plain; charset="iso-8859-1" Ettore, Difficult to 
> say why lsodes or any of the BDF solvers does not solve your problem, 
> as I don't know what your model looks like. However, my guess is that 
> your model is not a set of ordinary differential equations (ODEs) but 
> rather a differential algebraic equation problem (DAE). I suspect this 
> due to the matlab setting "massingular" = yes and "mass = M". If it is 
> a DAE, then you can use function radau from the R-package deSolve. 
> Function radau requires the DAE to be specified with a mass matrix M: 
> M y' = f(x, y) so I think this would be the integrator you need. An 
> alternative is to use R-function daspk. Hope this helps, Karline

	[[alternative HTML version deleted]]

