Solving ODEs - only positive solutions

Question

I am trying to to solve ODEs restricted to positive solutions, i.e.:

dx/dt=f(x)

with x>=0.

In MATLAB this is very easy to implement. Is there any workaround or package for R to restrict the solution-space to positive values only?

This is very crucial in my case and unfortunately there is no alternative. I searched for a while now but without any success. :-(

Answer 1

There's still not quite enough to go on here. For the sorts of problems I'm familiar with, modifying the system to operate on the scale of the log-transformed state variables works well (you can always back-transform the results e.g. to compare them with data). I have used this, for example, with the SIR model in epidemiology. I'm going to try with @MauritsEver's example, to illustrate transforming the system to operate on the log scale:

library(deSolve)
model <- function (time, y, parms) {
   with(as.list(c(y, parms)), {
       dlogN <-   r * (1 - exp(logN) / K)
       list(dlogN)
   })
}

# Starting conditions
y <- c(logN = log(0.1))
parms <- c(r = 0.1, K = 10)
times <- seq(0, 100, 1)
out <- as.data.frame(ode(y, times, model, parms))
out_backtran <- transform(out,N=exp(logN))
plot(N~time,data=out_backtran)

This approach has the following disadvantages:

• it won't handle solutions that are exactly on the boundary, and will have trouble with solutions that approach the boundary "too fast" (i.e. where state variables converge to zero in finite time)
• as written it requires manual translation. It would be entirely possible to write a system that allowed the user to enter a set of equations and a set of transformations and applied the transformations automatically, but it would be some effort.
• It may increase computational effort slightly (any time we have to use the original-scale value of the state variable we have to exponentiate)

Answer 2

Without any specific example code or details on the ODE it's difficult to be more specific. It could be quite simple, depending on the problem.

Here is a trivial example using deSolve and its function deSolve::subset.

# Example straight from the deSolve manual
library(deSolve);
model <- function (time, y, parms) {
    with(as.list(c(y, parms)), {
        dN <-   r * N * (1 - N / K);
        list(dN)
    })
}

# Starting conditions
y <- c(N = 0.1);
parms <- c(r = 0.1, K = 10);
times <- seq(0, 100, 1);

# Solve ODE and plot
out <- ode(y, times, model, parms);
plot(out, type = "l", xlim = c(0, 100));

We now impose a constraint on time and subset the solution.

# Constrain: time > 20 and plot
out.constrained <- subset(out, select = c("time", "N"), subset = time > 20);
plot(out.constrained, type = "l", xlim = c(0, 100));