% SIMIODESample.tex and simiode.cls for LaTeX2
% Version 11 August 2013
% 
% This class (simiode.cls) is a modification of that used by
% the journal PRIMUS - Problems, Resources, and Issues in
% Mathematics Undergraduate Studies.
%
% See http://www.tandfonline.com/loi/upri .
% 
% The class was developed by  Marc Frantz with suggestions of the 
% PRIMUS Editorial Board, especially John Thoo.
%
% This template is to be used with the  class file "simiode.cls". 
%  
% The latest version of this class and template files 
% can be found at
% 
%       <http://www.SIMIODE.org> in the Author section.
%
% A "Quick Start Guide" appears below.
%
% For information on using TeX or LaTeX, go to %                 <http://www.tug.org>.

% This file must be processed with version LaTeX2e or higher.
%
\documentclass{simiode}
\usepackage{amsmath,amssymb}
\usepackage{color}
\usepackage{graphicx}
%\usepackage{tikz,pgfplots}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Toggle the following line to turn the teacher version on (1) or off (0).
\def\TeacherEdition{0}
\newcommand{\teacher}[1]{\ifnum\TeacherEdition=1 #1 \fi}
% added by E. Sullivan:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% NOTE regarding double spacing %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Manuscripts should be submitted doubled spaced.  This is 
% accomplished via the \baselinestretch command below.
% However, you may obtain an idea of how your article will 
% look in final form, so that you can better plan your 
% presentation.  To see your manuscript formatted very 
% similarly to the on-line published format of  SIMIODE, 
% comment the following line out by inserting a percent sign % (%) at its beginning.

\renewcommand{\baselinestretch}{1.3}

%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% NOTE regarding macros %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Any macros defined for your paper should be contained in 
% the top matter. Likewise, any environment definitions such 
% as \newtheorem or \newenvironment should also go in the 
% top matter. 
% 
% Do not \input your macros as separate files within your 
% final version because more files make it harder for the 
% editors and users to keep track of and/or modify and 
% append your paper.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Top matter is above and 
% the beginning of the document below.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\begin{document}
\newpage
\title{\ifnum\TeacherEdition=1 TEACHER VERSION \else STUDENT VERSION \fi \\From Data to Differential Equations}



\markboth{From Data to Differential Equations}{From Data to Differential Equations}
% \markboth{Name of Scenario for Left Page running Heading}{Name of Scenario for Right Page running Heading}







 \author{Eric Sullivan and Kelly Cline\\
 Department of Mathematics, Computer Science, and Engineering \\ 
 Carroll College\\
 Helena MT 59625 USA\\
esullivan@carroll.edu  \, \& \, kcline@carroll.edu
}
\makeStitlePDFLaTex


\begin{abstract}
    Students  use their knowledge of second-order linear differential equations in
    conjunction with physical intuition of spring-mass systems to estimate the damping
    coefficient and spring constant from data.  The data is presented as {\it total
    distance traveled} instead of displacement so the students must rely on previous
    calculus knowledge as well as knowledge of differential equations and numerical
    differentiation.
\end{abstract}



\section*{SCENARIO DESCRIPTION}

A motion sensor was set up to measure the motion in a spring-mass system, but something
went wrong and the motion sensor measured the total distance traveled instead of simply
measuring the distance from the sensor.  The experiment was conducted three times (same
spring and same mass) with different initial conditions.  The total distance traveled is
given as the data sets in the Excel spreadsheet titled {\it SpringMassData.xlsx}.  The
conditions for the three experiments were:
\begin{description}
    \item[Data Set \#1:] Initial position: $y(0)=1$.  Initial velocity: $y'(0)=0$.
    \item[Data Set \#2:] Initial position: $y(0)=0.5$.  Initial velocity: $y'(0)=1$.
    \item[Data Set \#3:] Initial position: $y(0)=-0.75$.  Initial velocity: $y'(0)=-2.5$.
\end{description}

\noindent {\bf Experimental Setup:}
\begin{itemize}
    \item The sensor gathered data at a rate of 20 samples per second.
    \item The experiment was run for 5 seconds.
    \item $y(t)$ is the distance from the equilibrium position of the spring-mass system.
        Positive values of $y(t)$ indicate that the mass was above the equilibrium position.
        Negative values of $y(t)$ indicate that the mass was below the equilibrium position.
    \item There is some noise in the data.
\end{itemize}

\noindent {\bf Your Task:} Use the data in the associated Excel spreadsheet to determine
the governing differential equation along with estimates of the parameters. Include all
necessary plots and mathematics. 

\noindent {\bf Recall:} 
\begin{itemize}
    \item For an un-driven spring-mass system you can either use
        \[ m y''(t) + by'(t) + ky(t) = 0 \]
        where $m$ is the mass, $b$ is the damping coefficient, and $k$ is the spring
        constant or you can divide through by the mass and use
        \[ y''(t) + \delta y'(t) + \kappa y(t) = 0 \]
        where $\delta = b/m$ and $\kappa = k/m$.
    \item Recall that if $y(t)$ is displacement and $v(t) = y'(t)$ is velocity then the
        total distance traveled is given by the function
        \[ d(t) = \int_0^t \left| v(\tau) \right| d\tau \]
\end{itemize}


\teacher{
\section*{NOTES FOR TEACHERS}
This modeling scenario might work best as a team project in a differential equations
course or a numerical methods course.  The students will need to have a handle on
numerical differentiation at a minimum, as well as some basic principles of parameter
estimation.  What follows in this document is one approach to  parameter estimation.
The present approach is more geometric in nature and relies on a bit of educated guesswork
rather than more complex analytical tools.  Published together with this scenario is
a Mathematica  notebook that offers a more thorough parameter estimation.

Students will likely plot the data in a spreadsheet  first so we'll include the
plot here (see Figure \ref{fig:3-60-InitialDataPlot}). From these plots the students should deduce
that the spring-mass system in question comes to rest after a fairly short amount of time.
This should give clues to the student that the damping force is nontrivial.  
\begin{figure}[ht!]
    \begin{center}
        \includegraphics[width=0.5\columnwidth]{3-60-InitialDataPlot.eps}
    \end{center}
\vspace*{-.3 true cm}
    \caption{Plots of the original data}
    \label{fig:3-60-InitialDataPlot}
\vspace*{-.3 true cm}
\end{figure}


We'll assume that the differential equation is 
\begin{flalign}
    y'' + \delta y' + \kappa y = 0
    \label{eqn:governing_ode}
\end{flalign}
where $\delta = b/m$ is the damping coefficient divided by the mass of the object and
$\kappa = k/m$ is the spring constant divided by the mass of the object.  {\bf For the sake of
immediate satisfaction of any instructor, the data for this problem was built with $\delta =2$ and
$\kappa=10$.}

What follows is only one possible technique for finding the $\delta$ and $\kappa$.  The
technique demonstrated here is quite crude; leaving plenty of room for clever colleagues
and students to do better work.  The algebraic analysis, on the other hand, may take
significant prompting for many students.


The characteristic polynomial of  (\ref{eqn:governing_ode}) is $p(\lambda)=
\lambda^2 + \delta \lambda + \kappa$ and the roots of $p(\lambda)$ are given by the
quadratic formula as
\begin{flalign}
    \lambda_{1,2} = \frac{-\delta \pm \sqrt{\delta^2 - 4\kappa}}{2} =
    -\frac{\delta}{2} \pm \sqrt{\frac{\delta^2-4\kappa}{4}}
    \label{eqn:lambda12}
\end{flalign}

In the present situation it is clear that the spring-mass system eventually comes to rest.
Assuming that the damping and spring coefficients are such that the characteristic
polynomial has complex roots we can define $a=-\delta/2$ and $b = \sqrt{|\delta^2-4k|/4}$
so that $\lambda_{1,2} = a \pm b i$ and the solution to  (\ref{eqn:governing_ode})
is
\begin{flalign}
    y(t) = e^{at} \left( C_1 \sin(bt) + C_2 \cos(bt) \right).
    \label{eqn:gen_soln}
\end{flalign}

Knowing the basic setup for the solutions tells us two things:
\begin{enumerate}
    \item The values of $C_1$ and $C_2$ will depend on the initial conditions.
    \item The values of $a$ and $b$ will be the same for all runs of the experiment.
\end{enumerate}
Furthermore, we know that 
\[ y'(t) = b e^{at} \left( C_1\cos(bt) - C_2\sin(bt)  \right) + a e^{at} \left( C_1
    \sin(bt) + C_2 \cos(bt) \right) \]
so \vspace*{-.3 true cm}
\begin{subequations}
    \begin{flalign}
        y(0) &= C_2, \text{ and } \\
        y'(0) &= b C_1 + a C_2.
    \end{flalign}\label{eqn:constants}
\end{subequations}
From (\ref{eqn:constants}) we can write some fundamental relationships for each
of the experiments.
\begin{subequations}
    \begin{flalign}
        & \text{Experiment \#1: } C_2 = 1 \text{ and } C_1 = -\frac{a}{b} \quad \implies
        \quad y_1(t) = e^{at} \left( -\frac{a}{b} \sin(bt) + \cos(bt) \right)
        \label{eqn:constraints_specific_1} \\
        & \text{Experiment \#2: } C_2 = \frac{1}{2} \text{ and } C_1 = \frac{2-a}{2b}
        \quad \implies \quad y_2(t) = e^{at} \left( \frac{2-a}{2b} \sin(bt) + \frac{1}{2} \cos(bt) \right) \label{eqn:constraints_specific_2} \\
        & \text{Experiment \#3: } C_2 = -\frac{3}{4} \text{ and } C_1 = \frac{3a-10}{4b}
        \quad \implies \quad y_3(t) = e^{at} \left( \frac{3a-10}{4b} \sin(bt) - \frac{3}{4} \cos(bt) \right) \label{eqn:constraints_specific_3} 
    \end{flalign}
    \label{eqn:constraints_specific}
\end{subequations}

It should be clear that working with experiment \#1 and 
\eqref{eqn:constraints_specific_1} will likely be the simplest.  After the data from
experiment \#1 has been used to determine $a$ and $b$ we can check experiments \#2 and \#3
using \eqref{eqn:constraints_specific_2} and \eqref{eqn:constraints_specific_3}.  It
remains, then, to find approximations for $a$ and $b$ from experiment \#1.

Let's establish some notation.  Let $y(t)$ be the position of the mass at time $t$.
Therefore, $y'(t)$ is the velocity of the mass at time $t$ and the function
\begin{flalign}
    d(t) = \int_0^t |y'(\tau)| d\tau
    \label{eqn:total_dist_integral}
\end{flalign}
expresses the total distance traveled. The absolute value in 
    \eqref{eqn:total_dist_integral} will likely cause the students some trouble! 
Hence, the derivative of $d(t)$ will be the absolute value of the velocity: $d'(t) =
|y'(t)|$. Since we only have data we need to perform a numerical derivative of data set
\#1. This will result in a slight loss of information and the noise in the data doesn't
actually reveal an absolute value since some of the divided differences are negative.  The
numerical derivative of data set \#1 is shown in Figure 2.

\begin{figure}[ht!] \label{fig:3-60-NumDerivDataSet1}
    \begin{center}
        \includegraphics[width=0.5\columnwidth]{3-60-NumDerivDataSet1.eps}  
\vspace*{-.3 true cm}
    \caption{Numerical derivative of data set \#1.}
  \end{center}
\vspace*{-.3 true cm}
\end{figure}


At this point we will illustrate two possible ways to approximate the parameters $a$ and
$b$ from the data set: graphical way and using the Excel Solver to minimize residual errors.

\subsubsection*{Estimating Parameters Graphically}
Since this derivative comes from data set \#1 we can use the derivative of 
\eqref{eqn:constraints_specific_1},
\begin{flalign}
    y_1'(t) = - \left( \frac{b^2-a^2}{b} \right) e^{at} \sin(bt),
    \label{eqn:deriv_y1}
\end{flalign}
to make approximations of the period (or frequency) and the decay rate. 
The observant student may notice that the first {\it hump} changes direction at $t\approx1.05$;
hence giving an approximate period of $T\approx2.10$. This means that $b\approx 2\pi/2.10 \approx
2.992$.  Now it remains to approximate the rate of exponential decay: $a$. Given the
approximate period of approximately $2.10$ we can take the values at $t \approx 0.525,
1.575, 2.625, 3.675,$ and $4.725$ and perform an exponential regression on these points
to determine $a$. These points are chosen since they {\it should} correspond to
approximate maximum values of $|y'(t)|$. None of these $t$ values are in the data set so we'll have to settle
for approximations.  Furthermore, there is so much noise in the tail end of the data that if
we simply use the first three approximate times we should be {\it close} to the decay
rate. From the data for the numerical derivative we find the approximate maximum values
are $(0.525, 1.8)$, $(1.575, 0.70)$, and $(2.625, 0.30)$ with an exponential regression
(from Excel) of $m(t) \approx 2.7719 e^{-0.853 t}$.  Hence, $a \approx -0.853$.

Finally we can put all of the analysis together to determine the values of $\delta$ and
$\kappa$.
\begin{flalign}
    \notag a &= -\frac{\delta}{2} \approx -0.853 \\
    &\implies \delta \approx 1.706 \\
    %
    \notag b &= \sqrt{\frac{|\delta^2 - 4\kappa|}{4}} \\
    &\implies \kappa \approx 9.67963.
\end{flalign}
Therefore, from the data at hand, the differential equation is approximately
\[ y''(t) + 1.706 y'(t) + 9.67963 y(t) = 0. \]
The data for this problem was built with $\delta = 2$ and $\kappa = 10$ so these estimates
(however crude) are {\it close}.  



\subsubsection*{Estimating Parameters with Excel's Solver}
A second method for estimating the parameters is in the Excel file
\texttt{Data\_To\_ODE\_Excel\_Soln.xlsx}. The idea is that the numerical derivative of the
data, $d'(t)$, should roughly match $|y'(t)|$ from \eqref{eqn:deriv_y1}.  Using Excel we
calculate 
\begin{flalign}
    S = \sum_j \left( d'(t_j) - |y'(t_j)| \right)^2
    \label{eqn:ssr}
\end{flalign}
where $j$ counts over the data points.  Using Excel's Solver we minimize $S$ while varying
$a$ and $b$, and find that the minimum sum of the squares errors for Data Set \#1 is
approximately 3.96. Figure 3 shows a comparison between the numerical derivatives in
the data as along with a fitted form of $|y'(t)|$ where $a \approx -0.858$ and $b \approx
3.091$.
\begin{figure}[ht!]\label{fig:3-60-SolverImage}
    \begin{center}
        \includegraphics[width=0.5\columnwidth]{SolverImage.eps}
    
   \parbox{9 true cm}{ \caption{A numerical computation of $d'(t)$ along with a $|y'(t)|$ fitted using the 
    Excel's Solver ($a \approx -0.858$ and $b\approx 3.091$).}}
\end{center}

\vspace*{-.5 true cm}
\end{figure}


Plots of the data sets, along with the approximations, are shown in
Figures 4, 5, and
6.  The instructor should keep in mind that there are many
other ways to perform that data fit. Such analysis by more advanced students could
include more thorough parameter estimation with confidence intervals; perhaps leading to
upper and lower bounds on their estimates.


\begin{figure}[ht!]
    \begin{center}
        \includegraphics[width=0.5\columnwidth]{DataSet1Soln.eps}
 
    \parbox{9 true cm}{ \caption{Plot showing how the approximations compare to the function that was used to
    generate data set \#1.}}
   \end{center}
    \label{fig:3-60-DataSet1Soln}
\end{figure}

\begin{figure}[ht!]
    \begin{center}
        \includegraphics[width=0.5\columnwidth]{DataSet2Soln.eps}
    
    \parbox{9 true cm}{ \caption{Plot showing how the approximations compare to the function that was used to
    generate data set \#2.}}
\end{center}
    \label{fig:3-60-DataSet2Soln}
\end{figure}


\begin{figure}[ht!]
    \begin{center}
        \includegraphics[width=0.5\columnwidth]{DataSet3Soln.eps}
    \parbox{9 true cm}{ \caption{Plot showing how the approximations compare to the function that was used to
    generate data set \#3.}}
    \end{center}
    \label{fig:3-60-DataSet3Soln}
\end{figure}
}



\end{document}
