Comparison of ODE solvers

Alfonso R. Reyes

2017-11-09

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 <- [email protected]
        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 [email protected]@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