For the differential equation:
\[\dfrac{dy}{dt} = -5 \, e^{-t}\] the analytical solution is: \[y(t) = 5 \, e^{-t}\]
library(rODE)
# ODETest.R
setClass("ODETest", slots = c(
n = "numeric" # counts the number of getRate evaluations
),
contains = c("ODE")
)
setMethod("initialize", "ODETest", function(.Object, ...) {
.Object@n <- 0
.Object@state <- c(5.0, 0.0)
return(.Object)
})
setMethod("getExactSolution", "ODETest", function(object, t, ...) {
# analytical solution
return(5.0 * exp(-t))
})
setMethod("getState", "ODETest", function(object, ...) {
object@state
})
setMethod("getRate", "ODETest", function(object, state, ...) {
object@rate[1] <- -state[1]
object@rate[2] <- 1 # rate of change of time, dt/dt
object@n <- object@n + 1
object@rate
})
# constructor
ODETest <- function() {
odetest <- new("ODETest")
odetest
}## [1] "initialize"
## [1] "getExactSolution"
## [1] "getState"
## [1] "getRate"# This script can also be found under ./demo
# ComparisonRK45App.R
#
# Compares the solution by the RK45 ODE solver versus the analytical solution
ComparisonRK45App <- function(verbose = FALSE) {
ode <- new("ODETest")
ode_solver <- RK45(ode)
ode_solver <- setStepSize(ode_solver, 1)
ode_solver <- setTolerance(ode_solver, 1e-8)
time <- 0
while (time < 50) {
ode_solver <- step(ode_solver)
stepSize <- ode_solver@stepSize # update the step size
time <- time + stepSize
# ode <- ode_solver@ode
state <- getState(ode_solver@ode)
if (verbose)
cat(sprintf("time=%10f xl=%14e error=%14e n=%5d \n",
time, state[1],
(state[1] - getExactSolution(ode_solver@ode, time)),
ode_solver@ode@n))
}
cat("rate steps evaluated #", ode_solver@ode@n)
}
ComparisonRK45App(verbose = TRUE)## time= 0.063874 xl= 4.690617e+00 error= -4.288925e-11 n= 0
## time= 0.127748 xl= 4.400378e+00 error= -8.047341e-11 n= 0
## time= 0.191621 xl= 4.128097e+00 error= -1.132419e-10 n= 0
## time= 0.255495 xl= 3.872665e+00 error= -1.416458e-10 n= 0
## time= 0.319369 xl= 3.633037e+00 error= -1.661009e-10 n= 0
## time= 0.383243 xl= 3.408238e+00 error= -1.869878e-10 n= 0
## time= 0.447116 xl= 3.197347e+00 error= -2.046545e-10 n= 0
## time= 0.510990 xl= 2.999506e+00 error= -2.194187e-10 n= 0
## time= 0.574864 xl= 2.813907e+00 error= -2.315725e-10 n= 0
## time= 0.638738 xl= 2.639792e+00 error= -2.413807e-10 n= 0
## time= 0.702611 xl= 2.476451e+00 error= -2.490892e-10 n= 0
## time= 0.766485 xl= 2.323217e+00 error= -2.549192e-10 n= 0
## time= 0.830359 xl= 2.179464e+00 error= -2.590741e-10 n= 0
## time= 0.894233 xl= 2.044606e+00 error= -2.617395e-10 n= 0
## time= 0.958107 xl= 1.918093e+00 error= -2.630831e-10 n= 0
## time= 1.021980 xl= 1.799408e+00 error= -2.632583e-10 n= 0
## time= 1.085854 xl= 1.688067e+00 error= -2.624043e-10 n= 0
## time= 1.149728 xl= 1.583615e+00 error= -2.606482e-10 n= 0
## time= 1.213602 xl= 1.485626e+00 error= -2.581049e-10 n= 0
## time= 1.277475 xl= 1.393701e+00 error= -2.548779e-10 n= 0
## time= 1.341349 xl= 1.307463e+00 error= -2.510621e-10 n= 0
## time= 1.405223 xl= 1.226562e+00 error= -2.467428e-10 n= 0
## time= 1.469097 xl= 1.150666e+00 error= -2.419969e-10 n= 0
## time= 1.561338 xl= 1.079467e+00 error= 3.019212e-02 n= 0
## time= 1.653580 xl= 9.843494e-01 error= 2.753173e-02 n= 0
## time= 1.745822 xl= 8.976131e-01 error= 2.510576e-02 n= 0
## time= 1.838064 xl= 8.185195e-01 error= 2.289356e-02 n= 0
## time= 1.930306 xl= 7.463954e-01 error= 2.087628e-02 n= 0
## time= 2.022548 xl= 6.806264e-01 error= 1.903676e-02 n= 0
## time= 2.114789 xl= 6.206528e-01 error= 1.735933e-02 n= 0
## time= 2.207031 xl= 5.659637e-01 error= 1.582970e-02 n= 0
## time= 2.299273 xl= 5.160936e-01 error= 1.443486e-02 n= 0
## time= 2.391515 xl= 4.706178e-01 error= 1.316293e-02 n= 0
## time= 2.483757 xl= 4.291491e-01 error= 1.200307e-02 n= 0
## time= 2.575999 xl= 3.913345e-01 error= 1.094542e-02 n= 0
## time= 2.668240 xl= 3.568519e-01 error= 9.980957e-03 n= 0
## time= 2.760482 xl= 3.254077e-01 error= 9.101481e-03 n= 0
## time= 2.852724 xl= 2.967343e-01 error= 8.299501e-03 n= 0
## time= 2.944966 xl= 2.705874e-01 error= 7.568187e-03 n= 0
## time= 3.037208 xl= 2.467445e-01 error= 6.901313e-03 n= 0
## time= 3.129450 xl= 2.250025e-01 error= 6.293201e-03 n= 0
## time= 3.221692 xl= 2.051763e-01 error= 5.738673e-03 n= 0
## time= 3.313933 xl= 1.870971e-01 error= 5.233007e-03 n= 0
## time= 3.446050 xl= 1.706110e-01 error= 1.125464e-02 n= 0
## time= 3.578167 xl= 1.494959e-01 error= 9.861751e-03 n= 0
## time= 3.710284 xl= 1.309941e-01 error= 8.641246e-03 n= 0
## time= 3.842401 xl= 1.147821e-01 error= 7.571792e-03 n= 0
## time= 3.974518 xl= 1.005765e-01 error= 6.634696e-03 n= 0
## time= 4.106635 xl= 8.812897e-02 error= 5.813576e-03 n= 0
## time= 4.238752 xl= 7.722200e-02 error= 5.094079e-03 n= 0
## time= 4.370869 xl= 6.766489e-02 error= 4.463628e-03 n= 0
## time= 4.502986 xl= 5.929059e-02 error= 3.911203e-03 n= 0
## time= 4.635103 xl= 5.195269e-02 error= 3.427147e-03 n= 0
## time= 4.767220 xl= 4.552295e-02 error= 3.002998e-03 n= 0
## time= 4.899337 xl= 3.988896e-02 error= 2.631343e-03 n= 0
## time= 5.031454 xl= 3.495225e-02 error= 2.305684e-03 n= 0
## time= 5.163571 xl= 3.062650e-02 error= 2.020329e-03 n= 0
## time= 5.352264 xl= 2.683612e-02 error= 3.149051e-03 n= 0
## time= 5.540957 xl= 2.222140e-02 error= 2.607541e-03 n= 0
## time= 5.729650 xl= 1.840022e-02 error= 2.159150e-03 n= 0
## time= 5.918343 xl= 1.523612e-02 error= 1.787863e-03 n= 0
## time= 6.107037 xl= 1.261613e-02 error= 1.480423e-03 n= 0
## time= 6.295730 xl= 1.044667e-02 error= 1.225850e-03 n= 0
## time= 6.484423 xl= 8.650262e-03 error= 1.015054e-03 n= 0
## time= 6.673116 xl= 7.162767e-03 error= 8.405055e-04 n= 0
## time= 6.861809 xl= 5.931061e-03 error= 6.959726e-04 n= 0
## time= 7.050503 xl= 4.911159e-03 error= 5.762935e-04 n= 0
## time= 7.320555 xl= 4.066638e-03 error= 7.576659e-04 n= 0
## time= 7.590608 xl= 3.104224e-03 error= 5.783559e-04 n= 0
## time= 7.860661 xl= 2.369576e-03 error= 4.414816e-04 n= 0
## time= 8.130714 xl= 1.808790e-03 error= 3.370001e-04 n= 0
## time= 8.400767 xl= 1.380720e-03 error= 2.572453e-04 n= 0
## time= 8.670820 xl= 1.053958e-03 error= 1.963654e-04 n= 0
## time= 8.940872 xl= 8.045272e-04 error= 1.498934e-04 n= 0
## time= 9.210925 xl= 6.141271e-04 error= 1.144194e-04 n= 0
## time= 9.615977 xl= 4.687872e-04 error= 1.355111e-04 n= 0
## time= 10.021029 xl= 3.126539e-04 error= 9.037796e-05 n= 0
## time= 10.426081 xl= 2.085220e-04 error= 6.027679e-05 n= 0
## time= 10.831133 xl= 1.390720e-04 error= 4.020108e-05 n= 0
## time= 11.236185 xl= 9.275297e-05 error= 2.681176e-05 n= 0
## time= 11.817853 xl= 6.186084e-05 error= 2.500200e-05 n= 0
## time= 12.399521 xl= 3.457798e-05 error= 1.397517e-05 n= 0
## time= 12.981189 xl= 1.932785e-05 error= 7.811593e-06 n= 0
## time= 13.562857 xl= 1.080357e-05 error= 4.366385e-06 n= 0
## time= 14.439904 xl= 6.038807e-06 error= 3.360875e-06 n= 0
## time= 15.316950 xl= 2.512251e-06 error= 1.398206e-06 n= 0
## time= 16.193997 xl= 1.045141e-06 error= 5.816875e-07 n= 0
## time= 17.553692 xl= 4.347974e-07 error= 3.158107e-07 n= 0
## time= 18.913387 xl= 1.118798e-07 error= 8.133131e-08 n= 0
## time= 20.956872 xl= 2.878834e-08 error= 2.482998e-08 n= 0
## time= 23.000357 xl= 4.112791e-09 error= 3.599880e-09 n= 0
## time= 26.781657 xl= 5.875659e-10 error= 5.758751e-10 n= 0
## time= 30.562957 xl= 6.386756e-10 error= 6.384091e-10 n= 0
## time= 34.344257 xl= 6.942310e-10 error= 6.942249e-10 n= 0
## time= 38.125556 xl= 7.546190e-10 error= 7.546188e-10 n= 0
## time= 41.906856 xl= 8.202598e-10 error= 8.202598e-10 n= 0
## time= 45.688156 xl= 8.916104e-10 error= 8.916104e-10 n= 0
## time= 49.469456 xl= 9.691675e-10 error= 9.691675e-10 n= 0
## time= 53.250756 xl= 1.053471e-09 error= 1.053471e-09 n= 0
## rate steps evaluated # 0Notes. In this example, the number of iterations does not return from
ode_solver@ode@nthat is part of the classODETest. We will try to fix this in another example.
In this example, we create the environment object stack that will allow us to store temporary values or accumulators inside an S4 class.
library(rODE)
setClass("ODETest", slots = c(
stack = "environment" # environment object inside the class
),
contains = c("ODE")
)
setMethod("initialize", "ODETest", function(.Object, ...) {
.Object@stack$n <- 0 # "n" belongs to the class environment
.Object@state <- c(5.0, 0.0)
return(.Object)
})
setMethod("getExactSolution", "ODETest", function(object, t, ...) {
# analytical solution
return(5.0 * exp(-t))
})
setMethod("getState", "ODETest", function(object, ...) {
object@state
})
setMethod("getRate", "ODETest", function(object, state, ...) {
object@rate[1] <- -state[1]
object@rate[2] <- 1 # rate of change of time, dt/dt
object@stack$n <- object@stack$n + 1 # add 1 to the rate count
object@rate
})
# constructor
ODETest <- function() {
odetest <- new("ODETest")
odetest
}
# class implementation
ComparisonRK45App <- function(verbose = FALSE) {
ode <- new("ODETest")
ode_solver <- RK45(ode)
ode_solver <- setStepSize(ode_solver, 1)
ode_solver <- setTolerance(ode_solver, 1e-8)
cat(sprintf("%10s %14s %14s %5s \n", "time", "xl", "error", "n")) # header
time <- 0
while (time < 50) {
ode_solver <- step(ode_solver)
stepSize <- ode_solver@stepSize # update the step size
time <- time + stepSize
state <- getState(ode_solver@ode)
if (verbose)
cat(sprintf("%10f %14e %14e %5d \n",
time, state[1],
(state[1] - getExactSolution(ode_solver@ode, time)),
ode_solver@ode@stack$n))
}
cat("rate steps evaluated #", ode_solver@ode@stack$n)
}
ComparisonRK45App(verbose = TRUE)## [1] "initialize"
## [1] "getExactSolution"
## [1] "getState"
## [1] "getRate"
## time xl error n
## 0.063874 4.690617e+00 -4.288925e-11 16
## 0.127748 4.400378e+00 -8.047341e-11 22
## 0.191621 4.128097e+00 -1.132419e-10 28
## 0.255495 3.872665e+00 -1.416458e-10 34
## 0.319369 3.633037e+00 -1.661009e-10 40
## 0.383243 3.408238e+00 -1.869878e-10 46
## 0.447116 3.197347e+00 -2.046545e-10 52
## 0.510990 2.999506e+00 -2.194187e-10 58
## 0.574864 2.813907e+00 -2.315725e-10 64
## 0.638738 2.639792e+00 -2.413807e-10 70
## 0.702611 2.476451e+00 -2.490892e-10 76
## 0.766485 2.323217e+00 -2.549192e-10 82
## 0.830359 2.179464e+00 -2.590741e-10 88
## 0.894233 2.044606e+00 -2.617395e-10 94
## 0.958107 1.918093e+00 -2.630831e-10 100
## 1.021980 1.799408e+00 -2.632583e-10 106
## 1.085854 1.688067e+00 -2.624043e-10 112
## 1.149728 1.583615e+00 -2.606482e-10 118
## 1.213602 1.485626e+00 -2.581049e-10 124
## 1.277475 1.393701e+00 -2.548779e-10 130
## 1.341349 1.307463e+00 -2.510621e-10 136
## 1.405223 1.226562e+00 -2.467428e-10 142
## 1.469097 1.150666e+00 -2.419969e-10 148
## 1.561338 1.079467e+00 3.019212e-02 154
## 1.653580 9.843494e-01 2.753173e-02 160
## 1.745822 8.976131e-01 2.510576e-02 166
## 1.838064 8.185195e-01 2.289356e-02 172
## 1.930306 7.463954e-01 2.087628e-02 178
## 2.022548 6.806264e-01 1.903676e-02 184
## 2.114789 6.206528e-01 1.735933e-02 190
## 2.207031 5.659637e-01 1.582970e-02 196
## 2.299273 5.160936e-01 1.443486e-02 202
## 2.391515 4.706178e-01 1.316293e-02 208
## 2.483757 4.291491e-01 1.200307e-02 214
## 2.575999 3.913345e-01 1.094542e-02 220
## 2.668240 3.568519e-01 9.980957e-03 226
## 2.760482 3.254077e-01 9.101481e-03 232
## 2.852724 2.967343e-01 8.299501e-03 238
## 2.944966 2.705874e-01 7.568187e-03 244
## 3.037208 2.467445e-01 6.901313e-03 250
## 3.129450 2.250025e-01 6.293201e-03 256
## 3.221692 2.051763e-01 5.738673e-03 262
## 3.313933 1.870971e-01 5.233007e-03 268
## 3.446050 1.706110e-01 1.125464e-02 274
## 3.578167 1.494959e-01 9.861751e-03 280
## 3.710284 1.309941e-01 8.641246e-03 286
## 3.842401 1.147821e-01 7.571792e-03 292
## 3.974518 1.005765e-01 6.634696e-03 298
## 4.106635 8.812897e-02 5.813576e-03 304
## 4.238752 7.722200e-02 5.094079e-03 310
## 4.370869 6.766489e-02 4.463628e-03 316
## 4.502986 5.929059e-02 3.911203e-03 322
## 4.635103 5.195269e-02 3.427147e-03 328
## 4.767220 4.552295e-02 3.002998e-03 334
## 4.899337 3.988896e-02 2.631343e-03 340
## 5.031454 3.495225e-02 2.305684e-03 346
## 5.163571 3.062650e-02 2.020329e-03 352
## 5.352264 2.683612e-02 3.149051e-03 358
## 5.540957 2.222140e-02 2.607541e-03 364
## 5.729650 1.840022e-02 2.159150e-03 370
## 5.918343 1.523612e-02 1.787863e-03 376
## 6.107037 1.261613e-02 1.480423e-03 382
## 6.295730 1.044667e-02 1.225850e-03 388
## 6.484423 8.650262e-03 1.015054e-03 394
## 6.673116 7.162767e-03 8.405055e-04 400
## 6.861809 5.931061e-03 6.959726e-04 406
## 7.050503 4.911159e-03 5.762935e-04 412
## 7.320555 4.066638e-03 7.576659e-04 418
## 7.590608 3.104224e-03 5.783559e-04 424
## 7.860661 2.369576e-03 4.414816e-04 430
## 8.130714 1.808790e-03 3.370001e-04 436
## 8.400767 1.380720e-03 2.572453e-04 442
## 8.670820 1.053958e-03 1.963654e-04 448
## 8.940872 8.045272e-04 1.498934e-04 454
## 9.210925 6.141271e-04 1.144194e-04 460
## 9.615977 4.687872e-04 1.355111e-04 466
## 10.021029 3.126539e-04 9.037796e-05 472
## 10.426081 2.085220e-04 6.027679e-05 478
## 10.831133 1.390720e-04 4.020108e-05 484
## 11.236185 9.275297e-05 2.681176e-05 490
## 11.817853 6.186084e-05 2.500200e-05 496
## 12.399521 3.457798e-05 1.397517e-05 502
## 12.981189 1.932785e-05 7.811593e-06 508
## 13.562857 1.080357e-05 4.366385e-06 514
## 14.439904 6.038807e-06 3.360875e-06 520
## 15.316950 2.512251e-06 1.398206e-06 526
## 16.193997 1.045141e-06 5.816875e-07 532
## 17.553692 4.347974e-07 3.158107e-07 538
## 18.913387 1.118798e-07 8.133131e-08 544
## 20.956872 2.878834e-08 2.482998e-08 550
## 23.000357 4.112791e-09 3.599880e-09 556
## 26.781657 5.875659e-10 5.758751e-10 562
## 30.562957 6.386756e-10 6.384091e-10 568
## 34.344257 6.942310e-10 6.942249e-10 574
## 38.125556 7.546190e-10 7.546188e-10 580
## 41.906856 8.202598e-10 8.202598e-10 586
## 45.688156 8.916104e-10 8.916104e-10 592
## 49.469456 9.691675e-10 9.691675e-10 598
## 53.250756 1.053471e-09 1.053471e-09 604
## rate steps evaluated # 604