I have a problem integrating the 'standard map' (
http://en.wikipedia.org/wiki/Standard_maphttp://en.wikipedia.org/wiki/Standard_map ) with deSolve:
By using the modulo-operator '%%' with 2*pi in the ODEs (standardmap1), the
resulting values of P and Theta, should not be greater than 2pi. Because
this was not the case, i was thinking that the function 'ode' has some
internal problems with the '%%' or integrating periodical ODEs, so I wrote a
modulo-function by myself (modulo and standardmap2). But still I get values
much higher than 2pi and I cannot find the error... Any guess?
Thanks
Full code:
----------------------------------------------------------------
library(deSolve)
iterations <- 100
Parameter <- c(k = 0.6)
State <- c(Theta = 1 , P = 1)
Time <- 0:iterations/10
standardmap1 <- function(Time, State, Parameter){
with(as.list(c(State, Parameter)), {
dP <- (P + k * sin(Theta)) %% (2 * pi)
dTheta <- (P + Theta) %% (2 * pi)
return(list(c(dP, dTheta)))
})
}
#solve ode using standardmap1
out1 <- as.data.frame(ode(func = standardmap1, y = State,parms = Parameter,
times = Time))
# x mod y, end: maximal iterations
modulo <- function(x,y,end=1000){
for (n in 0:end)
if (x > (n-1)*y)
z = x - (n-1)*y
if (z > 0)
return(z)
else break
}
standardmap2 <- function(Time, State, Parameter){
with(as.list(c(State, Parameter)), {
dTheta <- modulo((P + Theta),(2*pi))
dP <- modulo(P + k *sin(Theta),(2*pi))
return(list(c(dP, dTheta)))
})
}
#solve ode using standardmap2
out2 <- as.data.frame(ode(func = standardmap2, y = State,parms = Parameter,
times = Time))
#plot the results
matplot(out1[2],out1[3], type = "p", pch = ".")
matplot(out2[2],out2[3], type = "p", pch = ".")
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-models at r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-dynamic-models
Hope it helps
Thomas Petzoldt