[R] deSolve - Function daspk on DAE system - Error (Vince)

[Soetaert, Karline]

Vince, 

When that happens, one possible reason is that your DAE is of index > 1,
which cannot be solved by daspk.

The solver radau, also from deSolve can handle DAEs up to index 3, but
you need to rewrite the problem in the form M*y' = f(x,y), where M is a
mass matrix.

If you do that for your problem, and solve it with radau, then radau
complains that the "matrix is repeatedly  singular", the problem is too
stiff, and it stops just like daspk. 
I think this means that this particular DAE is unsolvable, so you will
need to look at the formulation itself. 

By the way, there is a special R-mailinglist that deals with this type
of problems:
r-sig-dynamic-models at r-project.org


Hope this helps,


Karline

---------------------
Original message:

Date: Sat, 1 Oct 2011 20:20:10 -0700 (PDT)
From: Vince <vince.pileggi at ontario.ca>
To: r-help at r-project.org
Subject: [R] deSolve - Function daspk on DAE system - Error
Message-ID: <1317525610060-3864298.post at n4.nabble.com>
Content-Type: text/plain; charset=us-ascii

I'm getting this error on the attached code and breaking my head but
can't
figure it out. Any help is much appreciated. Thanks, Vince

CODE:
library(deSolve)

Res_DAE=function(t, y, dy, pars) {
  with(as.list(c(y, dy, pars)), {
     
  res1 = -dS -dES-k2*ES
  res2 = -dP + k2*ES
 
  eq1 = Eo-E -ES
  eq2 = So-S -ES -P
  return(list(c(res1, res2, eq1, eq2)))
  })
}

pars <- c(Eo=0.02, So=0.02, k2=250, E=0.01); pars
yini <- c(S=0.01, ES = 0.01, P=0.0, E=0.01); yini
times <- seq(0, 0.01, by = 0.0001); times
dyini = c(dS=0.0, dES=0.0, dP=0.0)

## Tabular output check of matrix output

DAE <- daspk(y = yini, dy = dyini, times = times, res = Res_DAE, parms =
pars, atol = 1e-10, rtol = 1e-10)

ERROR:
daspk--  warning.. At T(=R1) and stepsize H (=R2) the      nonlinear
solver
f
      nonlinear solver failed to converge      repeatedly of with abs
(H) =
H
      repeatedly of with abs (H) = HMIN      preconditioner had repeated
failur
        0.0000000000000D+00  0.5960464477539D-14

Warning messages:
1: In daspk(y = yini, dy = dyini, times = times, res = Res_DAE, parms =
pars,  :
  repeated convergence test failures on a step - inaccurate Jacobian or
preconditioner?
2: In daspk(y = yini, dy = dyini, times = times, res = Res_DAE, parms =
pars,  :
  Returning early. Results are accurate, as far as they go