Dear Albert2002,
there is no problem with deSolve and, of course, no problem with R's
modulo operator, but there are at least two errors in your model
formulation:
1.) The order of the returned derivatives must be exactly the same as
specified in the state variables. This is documented in the help files
and mentioned in section "Troubleshooting" (11.2) in the deSolve
vignette. You use:
State <- c(Theta = 1 , P = 1)
return(list(c(dP, dTheta)))
but you should use:
State <- c(P = 1, Theta = 1)
return(list(c(dP, dTheta)))
2.) Your model is **not a differential equation** but a difference
equation. It is (1) discrete (not continuous) and returns the new state
(not the derivative). As the name of deSolve suggests, this package is
primarily for differential equations. Nevertheless, it can be useful for
difference equations too, if one respects the distinction.
Solution A:
===========
Use the new development version of deSolve from
http://deSolve.r-forge.r-project.org
that has a solver method "iteration" for this type of models:
out1 <- ode(func = standardmap1, y = State,parms = Parameter,
times = Time, method = "iteration")
plot(out1)
Solution B:
===========
Rewrite your model so that it returns the 'derivative' and not the new
state and use method="euler". This works already with recent versions of
deSolve.
iterations <- 100
Parameter <- c(k = 0.6)
State <- c(P=1, Theta = 1 )
Time <- 0:iterations
standardmap2 <- function(Time, State, Parameter){
with(as.list(c(State, Parameter)), {
P1 <- (P + k * sin(Theta)) %% (2 * pi)
Theta1 <- (P + Theta) %% (2 * pi)
return(list(c(P1-P, Theta1-Theta)))
})
}
out2 <- ode(func = standardmap2, y = State, parms = Parameter,
times = Time, method = "euler")
plot(out2)
##------------------------------------------------------------
You may also consider to rewrite your problem in matrix form to get the
full map easily and without using loops.
If you have further questions, consider to subscribe to the "dynamic
models in R" mailing list:
r-sig-dynamic-mod...@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-dynamic-models
Hope it helps
Thomas Petzoldt
______________________________________________
R-help@r-project.org mailing list
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.
