[R] deSolve: Problem solving ODE including modulo-operator

[Thomas Petzoldt]

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