[R] Problem with ode

[Bingzhang Chen]

Hi,

I am trying to run a 1D nutrient-phytoplankton-zooplankton model in R
using the package 'deSolve'.  The code is shown below:

DEPTH = seq(2.5, 147.5, 5)
NPZ = function(t, state, params){
  with(as.list(params), {
     P <- state[1:NB]
     Z <- state[(NB + 1): (2*NB)]
     N <- state[(2*NB + 1): (3*NB)]
     F.I = function(z, hr){
                I0 = function(hr){
                        if (hr <=16){
                           450*sin(hr*pi/16)
                           } else {
                           0
                           }
                }
                 i = which(z == DEPTH)
                 I = I0(hr=hr)*0.43*exp(-z*Kw - Kc*(cumsum(delz*P)[i]
- delz*P[i]/2))
                 FI = alpha*I/(um^2 + (alpha*I)^2)^.5
          return(FI)
      }
        J1 = function(x){
        n = 10^3
        FI2 <- function(y)F.I(z = x, y)
        J1 = sum(sapply(runif(n, 0, 16), FI2))*16/24/n
       return(J1)
        }
        J = sapply(DEPTH, J1)         #the factor of light limitation
       deltaz <- c(rep(5, NB - 1), 2.5)
       P.Flux <- -D*diff(c(P,0))/deltaz
       Z.Flux <- -D*diff(c(Z,0))/deltaz
       N.Flux <- -D*diff(c(N,N0))/deltaz
       dP.dt = P*um*N/(N+Kn)*J - Z*Im*P^2/(P^2+Kp^2) - diff(c(0,P.Flux))/delz
       dZ.dt = Z*(eps*Im*P^2/(P^2+Kp^2) - g*Z) - diff(c(0,Z.Flux))/delz
       dN.dt = -P*um*N/(N+Kn)*J + Z*(1-eps)*Im*P^2/(P^2+Kp^2) + g*Z^2
-  diff(c(0,N.Flux))/delz
       return(list(c(dP.dt, dZ.dt, dN.dt)))
       })
}

params1 = c(um = 1,              #maximal phytoplankton growth rate, d-1
            Kn = 0.2,          #nitrogen half-saturation constant uM/L
            Kw = 0.04,         #light attenuation due to water, unit: m^-1
            Kc = 0.03,      #light attenuation by phytoplankton,
#unit: m^2 (mmol N)^-1
            alpha = 0.025,
            Im = 0.6,          #zooplankton maximal ingestion rate, unit: /d
            Kp = 1,            #zooplankton half-saturation for
ingestion, unit: mmol N/m^3
            eps = .3,          #Growth efficiency of zooplankton
            delz = 5,
            D = 2.6,             #diffusion coefficiency, unit: m^2/d
            g = 0.025,       #zooplankton mortality rate, unit: d^-1
(mmol N)^-1
            N0 = 14,         #Nitrate concentration at 150 m
            NB = length(DEPTH))        #Number of depth intervals
P.int = c(rep(0.07, 8),.12, .18,rep(0.21,3),rep(0.35,2),
rep(0.42,2),rep(0.07, 8), rep(0.01, 5))
Z.int = P.int/7
N.int = c(rep(0.1, 13), rep(4.55, 5), rep(10,12))
state = c(P.int, Z.int, N.int)
Time = seq(0,100, by = 1)
NPZ.out <- ode.1D(state, Time, NPZ, params1, nspec = 3, names = c("P","Z","N"))

After running for about 7 hours, the result returns:
DLSODA-  At current T (=R1), MXSTEP (=I1) steps
      taken on this call before reaching TOUT
In above message, I =
[1] 5000
In above message, R =
[1] 0.0498949
Warning：
1: In lsoda(y[ii], times, func = bmod, parms = parms, bandup = nspec *  :
  an excessive amount of work (> maxsteps ) was done, but integration
was not successful - increase maxsteps
2: In lsoda(y[ii], times, func = bmod, parms = parms, bandup = nspec *  :
  Returning early. Results are accurate, as far as they go

And the diagnostics(NPZ.out) shows:

lsoda return code
--------------------

  return code (idid) =  -1
  Excess work done on this call. (Perhaps wrong Jacobian type MF.)

--------------------
INTEGER values
--------------------

  1 The return code : -1
  2 The number of steps taken for the problem so far: 5000
  3 The number of function evaluations for the problem so far: 44438
  5 The method order last used (successfully): 1
  6 The order of the method to be attempted on the next step: 1
  7 If return flag =-4,-5: the largest component in error vector 0
  8 The length of the real work array actually required: 1732
  9 The length of the integer work array actually required: 110
 14 The number of Jacobian evaluations and LU decompositions so far: 4834
 15 The method indicator for the last succesful step,
           1=adams (nonstiff), 2= bdf (stiff): 2
 16 The current method indicator to be attempted on the next step,
           1=adams (nonstiff), 2= bdf (stiff): 2

--------------------
RSTATE values
--------------------

  1 The step size in t last used (successfully): 2.189586e-06
  2 The step size to be attempted on the next step: 1.039085e-05
  3 The current value of the independent variable which the solver has
reached: 0.0498949
  4 Tolerance scale factor > 1.0 computed when requesting too much accuracy: 0
  5 The value of t at the time of the last method switch, if any: 0.005314453


So I wonder what are the exact problems with this model. Is it just
because of inappropriate parameters? I also wonder how I can speed up
the calculation time in R.

Thanks,
-- 
Bingzhang Chen
Ph. D.,
State Key Lab of Marine Environmental Science,
College of Oceanography and Environmental Science,
Xiamen University,
Xiamen, Fujian 361005
P. R. China