In jmzobitz/demodelr: Simulating Differential Equations with Data

knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)

NOTE: THIS NEEDS TO BE UPDATED WITH THE REVISED CODE FROM 2021 Here we are going to show how we can utilize Euler's method to easily solve a system of differential equation. To recap, the most general form of a differential equation is:

$$ \displaystyle \frac{d\vec{y}}{dt} = f(\vec{y},\vec{\alpha},t), $$ where $\vec{y}$ is the vector of state variables you want to solve for, and $\vec{\alpha}$ is your vector of parameters.

At a given initial condition, Euler's method applies locally linear approximations to forecast the solution forward $\Delta t$ time units:

$$ \vec{y}{n+1} = y{n} + f(\vec{y}{n},\vec{\alpha},t{n}) \cdot \Delta t $$

Set up

In order to make Euler's method work, you will need four things:

• The size ($\Delta t$) of your timesteps.
• The number of timesteps you wish to run Euler's method. More timesteps means more computational time. If $N$ is the number of timesteps, $\Delta t \cdot N$ is the total time.
• A function that we have for our dynamics (this is $\displaystyle \frac{d\vec{y}}{dt} = f(\vec{y},\vec{\alpha},t)$),
• The values of the vector of parameters $\vec{\alpha}$.

We will work you through this step by step, with example code.

First you should open up a blank RStudio document. Be sure to have the library MAT369Code installed.

# Uncomment the next two lines if you need to either update or install the MAT 369 Code
# library(devtools)
# install_github("jmzobitz/MAT369Code",build_vignettes = TRUE)
library(MAT369Code)

Set up your time values

First you need to set up the length of each timestep ($\delta T$) and the number of timesteps we compute the solution. You do this by specifiying the variables deltaT and timeSteps:

deltaT <- 0.05    # timestep length
timeSteps <- 200   # must be a number greater than 1

Set up your parameters

Next you should specify the values of your parameters $\vec{alpha}$. The good news is that you can name the parameters and their values, and just refer to the shorthand name (not the numerical value) in subsequent code.

parameters <- c(r = 2, k = 0.5, e = 0.1, d = 1)   # parameters: a named vector

Set up your initial conditions.

Now we need to specify the starting values (initial conditions) for our differential equation. The code can easily accommodate multiple initial conditions (so you can test different initial conditions - I call this a "run").

initialCondition <- c(V=1, P=3)  # Be sure you have enough conditions as you do variables.

# To do multiple initial conditions, here is some sample code:
# initialCondition <- rbind(c(V=1, P=3),  # Separate out initial conditions with a separate row.
#                        c(V=5,P=6),
#                        c(V=10,P=9) )

Set up your dynamics

You need to have a function for the dynamics of the differential equation. The function euler looks for a function called dynamics.

The format of this function is structured, but easily modifiable. Here is a sample code that shows the dynamics for the Lotka-Volterra equations: $$ \begin{align} \frac{dV}{dt} &= r V - kVP \ \frac{dP}{dt} &= e k V P - dP \end{align} $$ For simplicity, I refer to the left hand side of the dynamics as dVariable so for example dV signifies $\displaystyle \frac{dV}{dt}$.

# R function to calculate the value of the derivatives at each time value
# Use the names of the variables as defined in the vectors above

dynamics <- function(t, state, parameters){
  with(as.list(c(state, parameters)), {  # Do not edit this line
    dV = r*V - k*V*P  # <-- You may edit this line
    dP = e*k*V*P - d*P # <-- You may edit this line.
    return(list(c(dV, dP)))  # <-- If you have more equations you will need to list the dVariables
  })
}

Putting it all together

Now we have everything we need! Basically the next step is to run the command euler, which will print to the screen a plot of your equation. If you have entered everything above correctly, you may copy the following code to your working file.

If you have multiple initial conditions, you will get a stacked plot for each run.

euler(deltaT,timeSteps,initialCondition,FUN=dynamics,parameters=parameters)

Your Turn!

• Try adjusting the number of solutions.
• Try putting in another differential equation and see what results you obtain.
• Try changing the step size.

jmzobitz/demodelr documentation built on March 6, 2024, 8:31 p.m.