#' DormandPrince45 class
#'
#' @param object a class object
#' @param enable a logical flag
#' @param ... additional parameters
#'
#' @rdname DormandPrince45-class
#'
#' @include ODEAdaptiveSolver.R ODE.R
#' @example ./inst/examples/KeplerDormandPrince45.R
setClass("DormandPrince45", slots = c(
error_code = "numeric",
a = "matrix",
b5 = "numeric",
er = "numeric",
numStages = "numeric",
stepSize = "numeric",
numEqn = "numeric",
temp_state = "numeric",
k = "matrix",
truncErr = "numeric",
ode = "ODE",
tol = "numeric",
enableExceptions = "logical"
),
contains = c("ODEAdaptiveSolver")
)
setMethod("initialize", "DormandPrince45", function(.Object, ode, ...) {
# initialized the Euler ODE solver
.Object@error_code <- .Object@NO_ERROR
.Object@a <- rbind( c(1.0/5.0, 0, 0, 0, 0),
c(3.0/40.0, 9.0/40.0, 0, 0, 0),
c(3.0/10.0, -9.0/10.0, 6.0/5.0, 0, 0),
c(226.0/729.0, -25.0/27.0, 880.0/729.0, 55.0/729.0, 0),
c(-181.0/270.0, 5.0/2.0, -266.0/297.0, -91.0/27.0, 189.0/55.0))
.Object@b5 <- c(19.0/216.0, 0.0, 1000.0/2079.0, -125.0/216.0, 81.0/88.0, 5.0/56.0)
.Object@er <- c(-11.0/360.0, 0.0, 10.0/63.0, -55.0/72.0, 27.0/40.0, -11.0/280.0)
.Object@numStages <- 6 # number of intermediate rate computations
.Object@stepSize <- 0.01
.Object@numEqn <- 0
.Object@ode <- ode # set the ode to ODESolver slot
.Object@tol <- 1.0e-6
.Object@enableExceptions <- FALSE
return(.Object)
})
#' @rdname init-method
setMethod("init", "DormandPrince45", function(object, stepSize, ...) {
# inititalize the solver
object@stepSize <- stepSize
state <- getState(object@ode)
if (is.null(state)) {
stop("state vector not defined")
return(object) # state vector not defined.
}
if (object@numEqn != length(state)) {
object@numEqn <- length(state)
object@temp_state <- vector("numeric", object@numEqn)
object@k <- matrix(data = 0, nrow = object@numStages, ncol = object@numEqn)
}
object@ode@state <- state
object
})
#' @rdname step-method
setMethod("step", "DormandPrince45", function(object, ...) {
object@error_code <- object@NO_ERROR
iterations <- 10
currentStep <- object@stepSize
error <- 0
state <- getState(object@ode)
object@k[1,] <- getRate(object@ode, state)
# NEW iteration
repeat {
iterations <- iterations - 1
currentStep <- object@stepSize
# compute the k's
cum <- 0
for (s in 2:object@numStages) {
for (i in 1:object@numEqn) {
object@temp_state[i] <- state[i]
for (j in 1:(s-1)) {
object@temp_state[i] <- object@temp_state[i] +
object@stepSize * object@a[s-1, j] * object@k[j, i]
cum <- cum + 1 # dummy cum to test how many iterations
}
}
# print k array
# get the ODE object, state and rate
# object@ode <- getRate(object@ode, object@temp_state, object@k[s,])
# assign rate to k vector
# object@k[s,] <- object@ode@rate # in Java rate is passed by param
object@k[s,] <- getRate(object@ode, object@temp_state) # in Java rate is passed by param
} # end for loop "s"
# compute the error
error <- 0
for (i in 1:object@numEqn) {
object@truncErr <- 0
for (s in 1:object@numStages) {
object@truncErr <- object@truncErr +
object@stepSize * object@er[s] * object@k[s, i]
}
error <- max(error, abs(object@truncErr))
}
if (error <= 1.4e-45) { # error too small to be meaningful,
error <- object@tol / 1.0e5 # increase step size x10
}
# find h step for the next try. Update stepSize
if (error > object@tol) { # shrink, no more than x10
fac <- 0.9 * (error / object@tol)^-0.25
object@stepSize <- object@stepSize * max(fac, 0.1)
} else if (error < object@tol / 10.0) { # grow, but no more than x10
fac <- 0.9 * (error / object@tol)^-0.2
if (fac > 1) { # sometimes fac is <1 because error/tol is close to 1
object@stepSize <- object@stepSize * min(fac, 10)
}
}
if (!((error > object@tol) && (iterations > 0))) break
} # end repeat loop
# advance the state
for (i in 1:object@numEqn) {
for (s in 1:object@numStages) {
state[i] <- state[i] + currentStep * object@b5[s] * object@k[s, i]
}
}
if (iterations == 0) {
object@error_code <- object@DID_NOT_CONVERGE
if (object@enableExceptions) {
warning("DormandPrince45 ODE SOlver did not converge")
}
}
object@ode@state <- state
return(object)
}
)
#' @rdname DormandPrince45-class
setMethod("enableRuntimeExceptions", "DormandPrince45", function(object, enable) {
object@enableExceptions <- enable
})
#' @rdname setStepSize-method
setMethod("setStepSize", "DormandPrince45", function(object, stepSize, ...) {
object@stepSize <- stepSize
object
})
#' @rdname getStepSize-method
setMethod("getStepSize", "DormandPrince45", function(object, ...) {
return(object@stepSize)
})
#' @rdname setTolerance-method
#' @example ./inst/examples/ComparisonRK45ODEApp.R
#' @family adaptive solver methods
setMethod("setTolerance", "DormandPrince45", function(object, tol) {
object@tol <- abs(tol)
if (object@tol < 1.0E-12) {
err_msg = "Error: Dormand-Prince ODE solver tolerance cannot be smaller than 1.0e-12." # $NON-NLS-1$
if (object@enableExceptions) {
stop(err_msg)
}
cat(err_msg, "\n")
object@tol <- 1.0e-12
}
return(object)
})
#' @rdname getTolerance-method
#' @family adaptive solver methods
setMethod("getTolerance", "DormandPrince45", function(object) {
return(object@tol)
})
#' @rdname getErrorCode-method
#' @family adaptive solver methods
setMethod("getErrorCode", "DormandPrince45", function(object) {
return(object@error_code)
})
#' DormandPrince45 constructor
#'
#' @param .ode an ODE object
#' @importFrom methods new
#' @export
#' @example ./inst/examples/KeplerDormandPrince45.R
DormandPrince45 <- function(.ode) {
dormandPrince45 <- new("DormandPrince45", .ode)
dormandPrince45 <- init(dormandPrince45, dormandPrince45@stepSize)
return(dormandPrince45)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.