In f0nzie/rODE: Ordinary Differential Equation (ODE) Solvers Written in R Using S4 Classes

In this vignette we will use the example from the book "Numerical Solutions of Ordinary Differential Equations" by Atkinson, Han and Stewart. Google Books

Build the ODE class (without class accumulator)

Comparison of solutions: RK45 vs analytical solution

For the differential equation:

$$\dfrac{dy}{dt} = - y$$ the analytical solution is: $$y(t) = e^{-t}$$

Find the errors if: $$y(0) = 1$$

library(rODE)

# ODETest.R

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(1.0, 0.0)
    return(.Object)
})

setMethod("getExactSolution", "ODETest", function(object, t, ...) {
    # analytical solution
    return(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
}

Euler

ComparisonEulerApp <- function(stepSize) {
    ode <- new("ODETest")
    ode_solver <- Euler(ode)
    ode_solver <- setStepSize(ode_solver, stepSize)
    time <-  0
    rowVector <- vector("list")
    i <-  1
    while (time < 5) {
        state <- getState(ode_solver@ode)
        time <- state[2]
        error <- getExactSolution(ode_solver@ode, time) - state[1]
        rowVector[[i]] <- list(step_size = stepSize, 
                               t = time, 
                               y_t = state[1], 
                               error = error, 
                               rel_err = error / getExactSolution(ode_solver@ode, time),
                               steps = ode_solver@ode@stack$n
                               )
        ode_solver <- step(ode_solver)
        i <- i + 1
    }
    data.table::rbindlist(rowVector)
}
dt <- ComparisonEulerApp(stepSize = 0.2)
dt[round(t,1) %in% c(1.0, 2.0, 3.0, 4.0, 5.0)]

# get a summary table for different step sizes
get_table <- function(stepSize) {
    dt <- ComparisonEulerApp(stepSize)
    dt[round(t,2) %in% c(1.0, 2.0, 3.0, 4.0, 5.0)]
}

step_sizes <- c(0.2, 0.1, 0.05)
dt_li <- lapply(step_sizes, get_table)
data.table::rbindlist(dt_li)

Using RK4

ComparisonRK4App <- function(stepSize) {
    ode <- new("ODETest")
    ode_solver <- RK4(ode)
    ode_solver <- setStepSize(ode_solver, stepSize)
    time <-  0
    rowVector <- vector("list")
    i <-  1
    while (time < 5) {
        state <- getState(ode_solver@ode)
        time <- state[2]
        error <- getExactSolution(ode_solver@ode, time) - state[1]
        rowVector[[i]] <- list(step_size = stepSize, 
                               t = time, 
                               y_t = state[1], 
                               error = error, 
                               rel_err = error / getExactSolution(ode_solver@ode, time),
                               steps = ode_solver@ode@stack$n
                               )
        ode_solver <- step(ode_solver)
        i <- i + 1
    }
    data.table::rbindlist(rowVector)
}

# get a summary table for different step sizes
get_table <- function(stepSize) {
    dt <- ComparisonRK4App(stepSize)
    dt[round(t,2) %in% c(1.0, 2.0, 3.0, 4.0, 5.0)]
}

step_sizes <- c(0.2, 0.1, 0.05)
dt_li <- lapply(step_sizes, get_table)
data.table::rbindlist(dt_li)

f0nzie/rODE documentation built on May 14, 2019, 10:34 a.m.