[R-sig-dyn-mod] deSolve::ode workaround for solutions with NaN

Stephen Paul Ellner spe2 at cornell.edu
Tue Jun 29 17:04:15 CEST 2010


The logs trick (as suggested by Woodrow Setzer) often takes care of things, 
in my experience. But it can fail if the state variable x(t)--> 0, so that
log x(t) eventually becomes -Inf and ode() then bombs. What I
have done in that case, is go back to non-log scale and intervene 
in the ODEs. If x(t)=0 is a boundary that can't be crossed without
trouble arising, the reason (usually) is that x(t)=0 implies dx/dt =0 but
ode() can overshoot the boundary and keep going. In your model that may
or may not be true, depending on parameter values. To prevent it,  
set the 'bad' entries of dQ (where overshoot has occurred) to zero:  
 
badQ <- (!is.finite(dQ))|(Q<=0) 
dQ[badQ] <- 0  

Here's an example: 

blowup=function(t,y,parms) {
	dy <- c( -abs(log(y[1])), -y[1]*y[2]); 
	if(parms){
		bady <- (!is.finite(dy))|(y<=0)
		dy[bady] <- 0 
	}
	return(list(dy))
}

out1 <- ode(y=c(.5,1),times=seq(0,2,by=0.1),func=blowup,parms=FALSE) # fails
out2 <- ode(y=c(.5,1),times=seq(0,2,by=0.1),func=blowup,parms=TRUE)  # runs 


Stephen P. Ellner (spe2 at cornell.edu)
Department of Ecology and Evolutionary Biology
Corson Hall, Cornell University, Ithaca NY 14853-2701
Phone (607) 254-4221    FAX (607) 255-8088

----------------------------------------------------------------------

Message: 1
Date: Mon, 28 Jun 2010 20:36:46 +0200
From: Maik Renner <maikrenner at googlemail.com>
To: r-sig-dynamic-models at r-project.org
Subject: [R-sig-dyn-mod] deSolve::ode workaround for solutions with
        NaN
Message-ID:
        <AANLkTin-QJbueRZe8fPa1NTHn6gQYWyiQrf8Vg7n9sIV at mail.gmail.com>
Content-Type: text/plain; charset=UTF-8

Dear List,

I have a question regarding unstable solutions using deSolve. Due to
my model equation I get NaN, when my state drops below zero.
Now, I wish to set the state to some default value, when this happens,
such that the simulation is continuing without NaNs.
However I could not figure out how to do that technically.

An if clause in the Model equation did not help. And I also tried the
event / root example in the vignette, but it did not run, having some
problems with my forcing inputs.

I have a model with one simple differential equation and one state and
a forcing vector. Here is the model implementation:

KirchSimeq19 <- function(t, state, parameters, input1, input2) {
 with(as.list(c(state, parameters)),{
  P <- input1(t)
  E <- input2(t)
  dQ <- exp(c1 + c2 * log(Q) + c3 * (log(Q))^2) * ( (P - k_E * E) / Q  - 1)
 list(c(dQ))
 }) # end with(as.list...
}



Any help  is appreciated,

Maik

--
Maik Renner, Dipl. Hydrologist
TU Dresden Germany
Tel.: +49 351 463-31341
http://tu-dresden.de/meteorologie



------------------------------

Message: 2
Date: Mon, 28 Jun 2010 14:51:06 -0400
From: Setzer.Woodrow at epamail.epa.gov
To: Special Interest Group for Dynamic Simulation Models in R
        <r-sig-dynamic-models at r-project.org>
Cc: r-sig-dynamic-models at r-project.org,
        r-sig-dynamic-models-bounces at r-project.org
Subject: Re: [R-sig-dyn-mod] deSolve::ode workaround for solutions
        with NaN
Message-ID:
        <OF77673D2F.272CA48C-ON85257750.00676D7D-85257750.00678E34 at epamail.epa.gov>

Content-Type: text/plain; charset=US-ASCII

What would happen if you re-wrote the system in terms of d log(Q) / dt =
(1/Q) dQ / dt?

R. Woodrow Setzer, Ph. D.
National Center for Computational Toxicology
http://www.epa.gov/comptox
US Environmental Protection Agency
Mail Drop B205-01/US EPA/RTP, NC 27711
Ph: (919) 541-0128    Fax: (919) 541-1194



More information about the R-sig-dynamic-models mailing list