Description Usage Arguments Examples
Defines an ODE object for any solver
ODE constructor
1 2 3 4 5 6 7 |
object |
a class object |
... |
additional parameters |
state |
current state |
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 | # ++++++++++++++++++++++++++++++++++++++++++++++++++ example: PendulumApp.R
# Simulation of a pendulum using the EulerRichardson ODE solver
suppressPackageStartupMessages(library(ggplot2))
importFromExamples("Pendulum.R") # source the class
PendulumApp <- function(verbose = FALSE) {
# initial values
theta <- 0.2
thetaDot <- 0
dt <- 0.1
pendulum <- Pendulum()
# pendulum@state[3] <- 0 # set time to zero, t = 0
pendulum <- setState(pendulum, theta, thetaDot)
pendulum <- setStepSize(pendulum, dt = dt) # using stepSize in RK4
pendulum@odeSolver <- setStepSize(pendulum@odeSolver, dt) # set new step size
rowvec <- vector("list")
i <- 1
while (getState(pendulum)[3] <= 40) {
rowvec[[i]] <- list(t = getState(pendulum)[3], # time
theta = getState(pendulum)[1], # angle
thetadot = getState(pendulum)[2]) # derivative of angle
pendulum <- step(pendulum)
i <- i + 1
}
DT <- data.table::rbindlist(rowvec)
return(DT)
}
# show solution
solution <- PendulumApp()
plot(solution)
# +++++++++++++++++++++++++++++++++++++++++++++++++ example: PendulumEulerApp.R
# Pendulum simulation with the Euler ODE solver
# Notice how Euler is not applicable in this case as it diverges very quickly
# even when it is using a very small `delta t``?ODE
importFromExamples("PendulumEuler.R") # source the class
PendulumEulerApp <- function(verbose = FALSE) {
# initial values
theta <- 0.2
thetaDot <- 0
dt <- 0.01
pendulum <- PendulumEuler()
pendulum@state[3] <- 0 # set time to zero, t = 0
pendulum <- setState(pendulum, theta, thetaDot)
stepSize <- dt
pendulum <- setStepSize(pendulum, stepSize)
pendulum@odeSolver <- setStepSize(pendulum@odeSolver, dt) # set new step size
rowvec <- vector("list")
i <- 1
while (getState(pendulum)[3] <= 50) {
rowvec[[i]] <- list(t = getState(pendulum)[3],
theta = getState(pendulum)[1],
thetaDot = getState(pendulum)[2])
pendulum <- step(pendulum)
i <- i + 1
}
DT <- data.table::rbindlist(rowvec)
return(DT)
}
solution <- PendulumEulerApp()
plot(solution)
# +++++++++++++++++++++++++++++++++++++++++++++++++++++++++ example KeplerApp.R
# KeplerApp solves an inverse-square law model (Kepler model) using an adaptive
# stepsize algorithm.
# Application showing two planet orbiting
# File in examples: KeplerApp.R
importFromExamples("Kepler.R") # source the class Kepler
KeplerApp <- function(verbose = FALSE) {
# set the orbit into a predefined state.
r <- c(2, 0) # orbit radius
v <- c(0, 0.25) # velocity
dt <- 0.1
planet <- Kepler(r, v) # make up an ODE object
solver <- RK45(planet)
rowVector <- vector("list")
i <- 1
while (getState(planet)[5] <= 10) {
rowVector[[i]] <- list(t = planet@state[5],
planet1.r = getState(planet)[1],
p1anet1.v = getState(planet)[2],
planet2.r = getState(planet)[3],
p1anet2.v = getState(planet)[4])
solver <- step(solver)
planet <- getODE(solver)
i <- i + 1
}
DT <- data.table::rbindlist(rowVector)
return(DT)
}
solution <- KeplerApp()
plot(solution)
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ base class: FallingParticleODE.R
# Class definition for application FallingParticleODEApp.R
setClass("FallingParticleODE", slots = c(
g = "numeric"
),
prototype = prototype(
g = 9.8
),
contains = c("ODE")
)
setMethod("initialize", "FallingParticleODE", function(.Object, ...) {
.Object@state <- vector("numeric", 3)
return(.Object)
})
setMethod("getState", "FallingParticleODE", function(object, ...) {
# Gets the state variables.
return(object@state)
})
setMethod("getRate", "FallingParticleODE", function(object, state, ...) {
# Gets the rate of change using the argument's state variables.
object@rate[1] <- state[2]
object@rate[2] <- - object@g
object@rate[3] <- 1
object@rate
})
# constructor
FallingParticleODE <- function(y, v) {
.FallingParticleODE <- new("FallingParticleODE")
.FallingParticleODE@state[1] <- y
.FallingParticleODE@state[2] <- v
.FallingParticleODE@state[3] <- 0
.FallingParticleODE
}
|
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.