[R-sig-dyn-mod] lsodes VS ode15s

[Soetaert, Karline]

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

________________________________

From: r-sig-dynamic-models-bounces at r-project.org on behalf of ettore.mosca at itb.cnr.it
Sent: Thu 29/12/2011 16:07
To: r-sig-dynamic-models at r-project.org
Subject: [R-sig-dyn-mod] lsodes VS ode15s



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', at 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 <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]]

_______________________________________________
R-sig-dynamic-models mailing list
R-sig-dynamic-models at r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-dynamic-models


-------------- next part --------------
A non-text attachment was scrubbed...
Name: not available
Type: application/ms-tnef
Size: 9241 bytes
Desc: not available
URL: <https://stat.ethz.ch/pipermail/r-sig-dynamic-models/attachments/20111229/5e063303/attachment.bin>