Soetaert, Karline
2011-Oct-03 07:18 UTC
[R] deSolve - Function daspk on DAE system - Error (Vince)
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
