Hi,

to diagnose your problem, I tried the model with your original parameters, using a Runge-Kutta fixed step integrator and smaller time steps. You may try to make them even smaller.

rk4 has no warranted accuracy, so it is less reliable than "lsoda" etc. However, it can be useful for debugging, because it runs constantly through, even if NA or NaN's occur. We see that the states grow drastically to extremely high (or even negative) values.

This points indeed to unrealistic parameter values, or to a model-misspecification, which means that necessary feedback mechanisms are missing, so that i12 can grow, even if nothing is left (i1, i2) from which it can grow from.

Regards, Thomas



init <- c(i1=10, i2=10, i12=0)
parameters <- c(alpha1=0.7, alpha2=0.5, beta1=0.5, beta2=0.3,
  gamma1=0.5, gamma2=0.5, delta=0.5, N=100)
times <- seq(0, 10, by=.01)

simul <- as.data.frame(ode(y = init, times = times, func = system,
  parms = parameters, method="rk4"))

head(simul, 200)
    time            i1            i2          i12
1   0.00  1.000000e+01  1.000000e+01 0.000000e+00
2   0.01  1.685572e+01  1.433028e+01 6.451953e-01
3   0.02  2.567008e+01  1.875164e+01 4.259949e+00
4   0.03  3.251911e+01  2.091655e+01 4.215409e+01
5   0.04 -1.270280e+00 -1.988876e+00 1.581789e+02
6   0.05 -8.108876e+02 -4.992572e+02 2.631133e+04
7   0.06 -4.166359e+73 -2.975971e+73 6.427467e+98
8   0.07           NaN           NaN          NaN
9   0.08           NaN           NaN          NaN

...



--
Dr. Thomas Petzoldt
Technische Universitaet Dresden
Faculty of Environmental Sciences
Institute of Hydrobiology
01062 Dresden, Germany

E-Mail: thomas.petzo...@tu-dresden.de
http://tu-dresden.de/Members/thomas.petzoldt

______________________________________________
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

Reply via email to