# Comparison of ODE solvers

## Build the ODE class (without class accumulator)

### Comparison of solutions: RK45 vs analytical solution

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"

## Build and run the application ComparisonRK45App

# 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 # 0

Notes. In this example, the number of iterations does not return from ode_solver@ode@n that is part of the class ODETest. We will try to fix this in another example.

## Storing the number of counts in a class environment object

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