We extend a project, suggested by Kurt Bryan in his new text Differential Equations: A Toolbox For Modeling The World, on solving numerically a variation of the Hill-Keller ODE:

\[v(t)' = P - kv(t)^r\]

where \(P\) is assumed to be a constant propulsive force per unit mass and \(v(t)\) is the velocity. This equation has no simple analytical solution for a fractional power (\(r=3/2\)) and it must be solved numerically. The project being presented goes beyond solving the ODE numerically and it is about fitting the model parameters \(P\) and \(k\) to Usain Bolt’s Olympic Data.

I implemented this project for our Mathematical Modeling II course, which I taught in spring 2021.

This course is designed for applied math majors in our applied math program. We typically have 25-30 students in the class, which is offered every spring, after the students take Mathematical Modeling I in the fall semester, and after taking numerical methods and a number of programming courses.

The biggest challenge was the computational component, given that it was a synchronous, fully online course. As I decided to use R, I used Zoom to get into my students’ computers and install for them R, LaTeX and RStudio on their personal computers. That way I made sure they all had the required software.

I also used the Rstudio Cloud, which offers a free option with limited resources and number of project hours per month. They do offer a discounted instructor’s account, which allows you to create a virtual classroom and share R Markdown projects in RStudio.

I made an effort to scaffold the projects and give the students examples with complete R code. I did use many examples from Kurt Bryan’s book, which I implemented in R Markdown and included them along with the complete R code in my lecture notes, which I also developed with R Markdown in RStudio.

A numerical solution can be obtained using either Euler’s method with a very small step or the improved Euler’s method or the Runge-Kutta’s method, which can be implemented from scratch using either R or Python in RStudio.

Using the pracma R package, we can visualize the direction field of the more general Hill-Keller model with \(r = 3/2\), \(P = 8.45\) and \(k = 0.205\) and initial condition \(v(0.165) = 0\). We also superimpose the Runge-Kutta solution using the rk4() solver from the pracma R package. .

On the next two slides we show the R code used to implement this visualization and the visualization itself.

library(pracma)
P<-8.45
r<-3/2
k<-0.205
f <- function(t,v) P - k*v^r
# time and velocity intervals
tt <- c(0.165,10) # in seconds
vv <- c(0,13) # in meters
# the direction field
vectorfield(f, tt, vv, scale=0.1, n=20, col="red")
sol <- rk4(f, a=tt[1], b=tt[2], y0=0, n=100)
lines(sol$x, sol$y, col="blue", lwd=2)

Next, we implement the linear interpolation \(v(t)\) of the discrete ODE solution approximation, obtained using the Runge-Kutta or Euler’s improvd method, for example.

For the linear interpolation, we use the connector() function from the mosaic R package, authored by Daniel Kaplan, Randall Pruim and Nicholas Horton. We also use the GitHub developmental version of the mosaicCalc R package for the plotting function slice_plot(). A continuous linear interpolation \(v(t)\) for a given dataset \((t_k,v_k)_{k=1}^n\) can be implemented by the R code below:

my_data <-data_frame(tk,vk)
v<-connector(vk ~ tk, data=my_data)

where vk represents the data vector \((v_k)_{k=1}^n\) and tk represents the data vector \((t_k)_{k=1}^n\), and my_data is an R dataframe with the two vectors vk and tk.

Once we have \(v(t)\) from the linear interpolation, we can integrate it with respect to time to find Bolt’s position \(x(t)\) as a function of time.

\[x(t,k,P) - x(t_0, k,P) =\int_{t_0}^t x'(s,k,P)ds = \int_{t_0}^t v(s,k,P)ds\]

Integrating the interpolated function v can be implemented in R using the code:

integral<-integrate(v, t0, t)
integral$value  # extract the integral value

Note that the initial conditions for this model come from Usain Bolt’s data, given on the next slide, and they are given by: \(t_0=0.165\), \(x(t_0)=0\) and \(v(t_0) = 0\).

These data are from Bolt’s Olympic Games’ 100 meter final in Beijing, 2008. Bolt’s times were a world record, which he lowered in 2009 to 9.58 seconds.

The data are in the form of (time, distance) pairs, where distance is measured in meters, horizontally along the track from the starting line, and time is measured in seconds elapsed from the firing of the starting gun.

The initial \((0.165, 0)\) data point indicates that Bolt had a reaction time of 0.165 seconds after the gun was fired before he started running and crossed the starting line.

Finally, we want fit the Hill-Keller two-parameter model for the sprinter’s distance \(x(t,P,k)\) to Bolt’s Olympic data by minimizing the sum of squared residuals \(S(k,P)\) as a function of the model parameters \(k\) and \(P\):

\[S(k,P) = \sum_{j=1}^{11} (x(t_j,k,P) - x_j)^2\] where \((t_j)_{j=1}^{11}\) is the vector of Bolt’s times and \((x_j)_{j=1}^{11}\) is the vector of Bolt’s positions for the given times.

The function \(S(k,P)\) quantifies the fit to the data given by any pair of parameters \(k\) and \(P\) in the model. The goal is to find that optimal pair \((k,P) = (k_0,P_0)\) that minimizes the function \(S(k,P)\).

# distance function from the general Hill-Keller's model
x <- function(time,k,P){
  # numerical solution with Runge-Kutta
  f <- function(t,v,new_P=P,new_k=k) new_P - new_k*v^(3/2)
  sol <- rk4(f, a=0.165, b=10, y0=0, n=100)
  vel_data <-data_frame(tk=sol$x,vk=sol$y)
  # linearly interpolated continuous velocity
  v<-connector(vk ~ tk, data=vel_data)
  # integral of velocity from t0 to time
  distance<-integrate(v, 0.165, time)
  return(distance$value)
}
x<-Vectorize(x)

To carry out the random search minimization we need good initial values, so it is helpful to start with a contour plot of \(S(k,P)\). However, we have to manually keep adjusting the ranges for \(k\) and \(P\) and keep zooming in on the surface until the global minimum is located and good initial approximations are found.

We can create a contour plot of the objective function S(k,P) using plotFun() from the mosaic package to help us find intervals for the initial approximations needed for the random search algorithm. The sum of squared residuals \(S(k,P) = \sum_{j=1}^{11} (x(t_j,k,P) - x_j)^2\) can be implemented in a vectorized way:

S <- function(k,P) sum((x(time=time,k,P) - distance)^2)

where time and distance are the vectors of Bolt’s time and position data, and x() is the distance function defined earlier.

# Bolt's time in seconds
time<-c(0.165, 1.85, 2.87, 3.78, 4.65, 5.50, 6.32, 7.14, 7.96, 8.79, 9.69)
# Bolt's position in meters
distance<-seq(0,100,by=10)
# Sum of squared residuals
S <- function(k,P) sum((x(time=time,k,P) - distance)^2)
S <- Vectorize(S) 
plotFun(S(k,P) ~ k & P, xlim=range(0.20,0.21), ylim=range(8.4,8.6))

We can plot the function as a surface in 3D. We can create an interactive 3D plot by using the interactive_plot() function from the mosaicCalc package.

The R code that implements the interactive surface plot is given below. Note that this is a wrapper function that creates a plotply graph whose interactivity is only available in HTML format that can be generated from the source R Markdown document in RStudio.

interactive_plot(S(k,P) ~ k & P, domain(k = range(0.20,0.210), P = range(8.4,8.5)))

It is hard to read values from a surface plot and the contour plots are much more useful for that purpose.

One approach to minimizing \(S(k,P)\) is to use multivariable calculus optimization by solving the equations:

\[\frac{\partial S}{\partial k} = 0, \quad \frac{\partial S}{\partial P} = 0\]

These two equations are nonlinear and rather complicated and to solve them we can use Newton’s method in two variables. However, it is easier to implement the random search algorithm to minimize a function of two variables. We can use small intervals around the rough estimates for the optimal values of \(k\) and \(P\) that we identified visually from the contour plot to implement the random search algorithm that minimizes the sum of squared residuals \(S(k,P)\).

The optimal values of the two parameters that minimize \(S(k,P)\) are:

\[k \approx 0.206, P \approx 8.467 \, \to S\approx 1.161\]

a<-0.20 # lower limit of k
b<-0.21 # upper limit of k
c<-8.4 # lower limit of P
d<-8.5 # upper limit of P
N<-1e4 # number of iterations
kmin<-runif(1,a,b) # random number from U(a,b)
Pmin<-runif(1,c,d) # random number from U(c,d)
Smin<-S(kmin,Pmin) # objective function to minimize
for (n in 1:N){
  k<-runif(1,a,b) # update k with a random number from U(a,b)
  P<-runif(1,c,d) # update P with a random number from U(c,d)
  Sval<-S(k,P) # update objective function
  if (Sval < Smin) {
    kmin<-k
    Pmin<-P
    Smin<-Sval
  }}

k<-0.205 # optimal value
P<-8.45 # optimal value
# Bolt's data and fitted Hill-Keller's model
gf_point(distance~time, size=2, col="blue", 
        xlab="time (seconds)", ylab="position (meters)") %>%
  slice_plot(x(time=t, k=k, P=P) ~ t, domain(t=range(0,10)), col="red")

This project presents a general approach to solving numerically an ODE for velocity, then interpolating the discrete numerical solution to obtain a continuous velocity function, which can be integrated numerically to obtain a continuous distance function that can be fitted to time and distance data by minimizing the sum of squared residuals and using a random search algorithm to find the optimal parameter values.

It is often the case that fitting models with more parameters may lead to only minor improvements to the model’s agreement with the data. In general, our goal is to find a model with the fewest number of parameters, which provides a reasonable fit to the data, without overfitting.

We implement everything from scratch in R, using dynamic R Markdown documents in RStudio, which allow us to unify plain text narrative, mathematical typesetting with LaTeX and computing in R and Python.