{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#A Sage script  to fit the Hill-Keller model v'(t) = P - k*v(t) to the data from\n",
    "#Usain Bolt's 2008 Beijing gold medal Olympic race.  This notebook fits both k and P\n",
    "#to fit the data, in a least-squares sense.\n",
    "\n",
    "#The Data: Here is the data for Usain Bolt's 2008 Beijing race. The data consists\n",
    "#of pairs [t, d], where d is 0, 10, 20, ..., 100 meters and t is the corresponding\n",
    "#split time.\n",
    "\n",
    "bolttimes = [0.165, 1.85, 2.87, 3.78, 4.65, 5.50, 6.32, 7.14, 7.96, 8.79, 9.69];\n",
    "dists = [10*i for i in range(11)];\n",
    "data = list(zip(bolttimes,dists));"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#A plot\n",
    "plt1 = scatter_plot(data,rgbcolor=[1,0,0]);\n",
    "show(plt1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Solving The ODE: Solve the Hill-Keller ODE, using P = 11 and treating \"k\" as\n",
    "#unknown. Also use initial condition v(0.165) = 0.\n",
    "var('k P t');\n",
    "v = function('v')(t)\n",
    "t0 = 0.165;\n",
    "de = diff(v,t) == P - k*v\n",
    "vsol(t) = desolve(de,v,[t0,0],t);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Integrate to obtain the position in terms of t, k, and P:\n",
    "var('tau');\n",
    "X(t,k,P) = integral(vsol(tau),tau,t0,t);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Let's guess a value of k and plot X(t) with the data.\n",
    "plt2 = plot(X(t,1,11), t, 0, 9.69);\n",
    "pp = plt1 + plt2;\n",
    "show(plt1+plt2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#The Optimal Choice for k and P: Not bad, but we can do better by forming a sum of squares SS and minimizing with respect to k and P.\n",
    "SS = function('SS')(k,P);\n",
    "SS(k,P) = add((X(bolttimes[i])-dists[i])^2 for i in range(11));"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#A quick plot of log(SS(k,P))\n",
    "plot3d(log(SS), (k,0.7,1.1), (P,10,11)).show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Somewhere around k = 0.9, P=10.5 it appears there is a minimum.\n",
    "#We could minimize using Calculus 3 techniques, but Sage has a\n",
    "#built-in command \"find_fit\" to find the least-squares model fit to the data:\n",
    "model(t) = X(t,k,P); #Specify the model to be fit, with \"t\" as the independent variable\n",
    "sol = find_fit(data,model,parameters=[k,P],initial_guess=[0.9,10.5]) #Fit the model by adjusting k and b\n",
    "kbest = sol[0].rhs();\n",
    "Pbest = sol[1].rhs();\n",
    "bestX(t) = model(k=kbest,P=Pbest) #Define the best fit function of the form in \"model\"\n",
    "print(kbest,Pbest)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Plot the function X(t) with this best choice for k and P.\n",
    "plt2 = plot(bestX(t), t, 0, 9.69);\n",
    "pp = plt1 + plt2;\n",
    "show(plt1+plt2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "SageMath 9.2",
   "language": "sage",
   "name": "sagemath"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}