R help - Jan 2011 - deSolve: Problem solving ODE including modulo-operator

I have a problem integrating the 'standard map' (
http://en.wikipedia.org/wiki/Standard_map
http://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 = ".")
-- 
View this message in context:
http://r.789695.n4.nabble.com/deSolve-Problem-solving-ODE-including-modulo-operator-tp3236396p3236396.html
Sent from the R help mailing list archive at Nabble.com.

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