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

The FallingParticleODE class using the EUler ODE solver

Because this is free fall of a particle in one-dimension, for vertical motion, we use Newton's second law:

$$ m \frac{d^2y} {dt^2} = F $$ and we know that the gravitational force is: $$ F = -mg $$ Therefore: $$ m \frac{d^2y} {dt^2} = -g $$ That expressing as first-order differential equations: $$ \frac {dy}{dt} = v \ \frac {dv}{dt} = -g $$

($y$), we define the numerical state equations as:

state[1]  x
state[2]  v
state[3]  t

From the equations of motion for a falling particle, the derivatives are: $$ \dot s_1 = s_2 \ \dot s_2 = -g \ \dot s_3 = 1 $$ which is equivalent of writing this as the rate in the code: $$ r_1 = r_2 \ r_2 = -g \ r_3 = 1 $$

# This code can also be found in the `examples` folder under this name:
# FallingParticleODE.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
}


# This code can also be found in the `examples` folder under this name:
# 
# FallingParticleODEApp.R
#
#
FallingParticleODEApp <- function(verbose = FALSE) {
    library(ggplot2)

    # load the R class that sets up the solver for this application

    initial_y <- 10   # initial y position
    initial_v <- 0    # initial x position
    dt        <- 0.01 # delta time for step


    ball <- FallingParticleODE(initial_y, initial_v)

    solver <- Euler(ball)
    solver <- setStepSize(solver, dt)

    rowVector <- vector("list")
    i <- 1
    # stop loop when the ball hits the ground
    while (ball@state[1] >= 0) {
        rowVector[[i]] <- list(state1 = ball@state[1], 
                               state2 = ball@state[2], 
                               state3 = ball@state[3])
        solver <- step(solver)
        ball <- solver@ode
        if (verbose) {
            cat(sprintf("%12f %12f ",  ball@state[1], ball@rate[1] ))
            cat(sprintf("%12f %12f ",  ball@state[2], ball@rate[2] ))
            cat(sprintf("%12f %12f\n", ball@state[3], ball@rate[3] ))
        }
        i <- i + 1
    }

    DT <- data.table::rbindlist(rowVector)
    print(ggplot(DT, aes(x = state3, y = state1)) + geom_line(col = "blue"))
    print(ggplot(DT, aes(x = state3, y = state2)) + geom_line(col = "red"))
}


FallingParticleODEApp()

Plots

The figure sizes have been customised so that you can easily put two images side-by-side.

# ggplot(DT, aes(x = state3, y = state1)) + geom_line(col = "blue")
# ggplot(DT, aes(x = state3, y = state2)) + geom_line(col = "red")

Data

You can write math expressions, e.g. $Y = X\beta + \epsilon$, footnotes^[A footnote here.], and tables, e.g. using knitr::kable().

# knitr::kable(head(DT, 10))

The FallingParticleODE application

The complete R script can be found in the examples folder of this package. At the bottom of the document

AlfonsoRReyes/rODE documentation built on May 20, 2019, 6:48 p.m.