Description Usage Arguments Examples
DormandPrince45 ODE solver class
DormandPrince45 generic
DormandPrince45 constructor ODE
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 | DormandPrince45(ode, ...)
## S4 method for signature 'DormandPrince45'
init(object, stepSize, ...)
## S4 replacement method for signature 'DormandPrince45'
init(object, ...) <- value
## S4 method for signature 'DormandPrince45'
step(object, ...)
## S4 method for signature 'DormandPrince45'
enableRuntimeExceptions(object, enable)
## S4 method for signature 'DormandPrince45'
setStepSize(object, stepSize, ...)
## S4 method for signature 'DormandPrince45'
getStepSize(object, ...)
## S4 method for signature 'DormandPrince45'
setTolerance(object, tol)
## S4 replacement method for signature 'DormandPrince45'
setTolerance(object, ...) <- value
## S4 method for signature 'DormandPrince45'
getTolerance(object)
## S4 method for signature 'DormandPrince45'
getErrorCode(object)
## S4 method for signature 'ODE'
DormandPrince45(ode, ...)
|
ode |
ODE object |
... |
additional parameters |
object |
a class object |
stepSize |
size of the step |
value |
step size to set |
enable |
a logical flag |
tol |
tolerance |
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 | # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ base class: KeplerVerlet.R
setClass("KeplerDormandPrince45", slots = c(
GM = "numeric",
odeSolver = "DormandPrince45",
counter = "numeric"
),
contains = c("ODE")
)
setMethod("initialize", "KeplerDormandPrince45", function(.Object, ...) {
.Object@GM <- 4 * pi * pi # gravitation constant times combined mass
.Object@state <- vector("numeric", 5) # x, vx, y, vy, t
.Object@odeSolver <- DormandPrince45(.Object)
.Object@counter <- 0
return(.Object)
})
setMethod("doStep", "KeplerDormandPrince45", function(object, ...) {
object@odeSolver <- step(object@odeSolver)
object@state <- object@odeSolver@ode@state
object
})
setMethod("getTime", "KeplerDormandPrince45", function(object, ...) {
return(object@state[5])
})
setMethod("getEnergy", "KeplerDormandPrince45", function(object, ...) {
ke <- 0.5 * (object@state[2] * object@state[2] +
object@state[4] * object@state[4])
pe <- -object@GM / sqrt(object@state[1] * object@state[1] +
object@state[3] * object@state[3])
return(pe+ke)
})
setMethod("init", "KeplerDormandPrince45", function(object, initState, ...) {
object@state <- initState
# call init in AbstractODESolver
object@odeSolver <- init(object@odeSolver, getStepSize(object@odeSolver))
object@counter <- 0
object
})
setReplaceMethod("init", "KeplerDormandPrince45", function(object, ..., value) {
object@state <- value
# call init in AbstractODESolver
object@odeSolver <- init(object@odeSolver, getStepSize(object@odeSolver))
object@counter <- 0
object
})
setMethod("getRate", "KeplerDormandPrince45", function(object, state, ...) {
# Computes the rate using the given state.
r2 <- state[1] * state[1] + state[3] * state[3] # distance squared
r3 <- r2 * sqrt(r2) # distance cubed
object@rate[1] <- state[2]
object@rate[2] <- (- object@GM * state[1]) / r3
object@rate[3] <- state[4]
object@rate[4] <- (- object@GM * state[3]) / r3
object@rate[5] <- 1 # time derivative
object@counter <- object@counter + 1
object@rate
})
setMethod("getState", "KeplerDormandPrince45", function(object, ...) {
# Gets the state variables.
return(object@state)
})
setReplaceMethod("setSolver", "KeplerDormandPrince45", function(object, value) {
object@odeSolver <- value
object
})
# constructor
KeplerDormandPrince45 <- function() {
kepler <- new("KeplerDormandPrince45")
return(kepler)
}
# +++++++++++++++++++++++++++++++++++++++++ Example: ComparisonRK45ODEApp.R
# Updates the ODE state instead of using the internal state in the ODE solver
# Also plots the solver solution versus the analytical solution at a
# tolerance of 1e-6
# Example file: ComparisonRK45ODEApp.R
# ODE Solver: Runge-Kutta 45
# ODE class : RK45
# Base class: ODETest
library(ggplot2)
library(dplyr)
library(tidyr)
importFromExamples("ODETest.R")
ComparisonRK45ODEApp <- function(verbose = FALSE) {
ode <- new("ODETest") # new ODE instance
ode_solver <- RK45(ode) # select ODE solver
ode_solver <- setStepSize(ode_solver, 1) # set the step
# two ways to set tolerance
# ode_solver <- setTolerance(ode_solver, 1e-6)
setTolerance(ode_solver) <- 1e-6
time <- 0
rowVector <- vector("list") # row vector
i <- 1 # counter
while (time < 50) {
# add solution objects to a row vector
rowVector[[i]] <- list(t = getState(ode)[2],
ODE = getState(ode)[1],
s2 = getState(ode)[2],
exact = getExactSolution(ode, time),
rate.counts = getRateCounts(ode),
time = time )
ode_solver <- step(ode_solver) # advance solver one step
stepSize <- getStepSize(ode_solver) # get the current step size
time <- time + stepSize
ode <- getODE(ode_solver) # get updated ODE object
state <- getState(ode) # get the `state` vector
i <- i + 1 # add a row vector
}
DT <- data.table::rbindlist(rowVector) # create data table
return(DT)
}
solution <- ComparisonRK45ODEApp()
plot(solution)
# aditional plot for analytics solution vs. RK45 solver
solution.multi <- solution %>%
select(t, ODE, exact)
plot(solution.multi) # 3x3 plot
# plot comparative curves analytical vs ODE solver
solution.2x1 <- solution.multi %>%
gather(key, value, -t) # make a table of 3 variables. key: ODE/exact
g <- ggplot(solution.2x1, mapping = aes(x = t, y = value, color = key))
g <- g + geom_line(size = 1) +
labs(title = "ODE vs Exact solution",
subtitle = "tolerance = 1E-6")
print(g)
|
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.