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

[Maik Renner]

Dear all,



thank you a lot for your help and suggestions.

I tried Stephens approach with some success and it was the solution I
was looking for.

I still have to play around with the values I set for dy.

Because setting it to zero, results in zeros for the remaining
simulation. But this is problem specific.

Maik
--

On Tue, Jun 29, 2010 at 5:04 PM, Stephen Paul Ellner <spe2 at cornell.edu> wrote:
> 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
>
> _______________________________________________
> R-sig-dynamic-models mailing list
> R-sig-dynamic-models at r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-dynamic-models
>



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