demo/KeplerDormandPrince45App.R In f0nzie/rODE: Ordinary Differential Equation (ODE) Solvers Written in R Using S4 Classes

```# ++++++++++++++++++++++++++++++++++++++++++ example: KeplerDormandPrince45App.R
# Demostration of the use of ODE solver RK45 for a particle subjected to
# a inverse-law force. The difference with the example KeplerApp is we are
# seeing the effect in thex and y axis on the particle.
# The original routine used the Verlet ODE solver

importFromExamples("KeplerDormandPrince45.R")

set_solver <- function(ode_object, solver) {
slot(ode_object, "odeSolver") <- solver
ode_object
}

KeplerDormandPrince45App <- function(verbose = FALSE) {
# values for the examples
x  <- 1
vx <- 0
y  <- 0
vy <- 2 * pi
dt <- 0.01          # step size
tol <- 1e-3         # tolerance
particle  <- KeplerDormandPrince45()                      # use class Kepler

# Two ways of initializing the ODE object
# particle  <- init(particle, c(x, vx, y, vy, 0))  # enter state vector
init(particle) <- c(x, vx, y, vy, 0)

odeSolver <- DormandPrince45(particle)      # select the ODE solver

# Two ways of initializing the solver
# odeSolver <- init(odeSolver, dt)            # start the solver
init(odeSolver) <-  dt

# Two ways of setting the tolerance
# odeSolver <- setTolerance(odeSolver, tol) # this works for adaptive solvers
setTolerance(odeSolver) <- tol
setSolver(particle) <-  odeSolver

initialEnergy <- getEnergy(particle)        # calculate the energy
rowVector <- vector("list")
i <- 1
while (getTime(particle) < 1.5) {
rowVector[[i]] <- list(t  = getState(particle)[5],
x  = getState(particle)[1],
vx = getState(particle)[2],
y  = getState(particle)[3],
vx = getState(particle)[4],
energy = getEnergy(particle) )
particle <- doStep(particle)            # advance one step
energy   <- getEnergy(particle)         # calculate energy
i <- i + 1
}
DT <- data.table::rbindlist(rowVector)
return(DT)
}

solution <- KeplerDormandPrince45App()
plot(solution)
```
f0nzie/rODE documentation built on May 14, 2019, 10:34 a.m.