Let $x(t)$ measure the distance the mass of a damped harmonic oscillator is from its equilibrium position at time $t$.  Note. The equilibrium position is assumed to be at $x=0$.  See image below.

The damped harmonic oscillator can be modeled by the following Differential Equation (DE) which we present as an Initial Value Problem (IVP) since we are including Initial Conditions (IC):

$m\ddot{x}+c\dot{x}+kx=0$
$x(0)=x_0,\dot{x}(0)=v_0$

$m=$  mass, $c=$  the damping constant, and $k=$ the spring constant, with $m,c,k>0$.

We divide the DE by $m$ to get the equivalent DE:

$\ddot{x}+\frac{c}{m}\dot{x}+\frac{k}{m}x=0$
$x(0)=x_0,\dot{x}(0)=v_0$

which is easier to work with. If we want, we can think of $c$ as representing $\frac{c}{m}$ and $k$ as representing $\frac{k}{m}$. With these substitutions the DE becomes:

$\ddot{x}+c\dot{x}+kx=0$
$x(0)=x_0,\dot{x}(0)=v_0$

The letters $c$ and $k$ appearing in the differential equation
$$\ddot{x}+c\dot{x}+kx=0$$
are called parameters.  Sometimes the values of parameters appearing in a DE can be measured directly. For example, if we go back to our first DE for the damped harmonic oscillator:
$$m\ddot{x}+c\dot{x}+kx=0$$
the parameter $m$, which is the mass, could be found by putting the mass on a scale and weighing it. However, often it is the case that there is no easy way to directly measure parameters. Instead their values need to be estimated from data, the data being measurements made on the entire system.

For example, with the damped harmonic oscillator we might have measured the position of the mass at several different times. This would result in our having a set of data points. One data point $(t,x)$ for each (time,  position) measurement:

$$(t_i,x_i)\text{ with }i=0,1,2,\dots ,n$$

This project. In this project we will estimate the parameters $c$ and $k$ of the DE

$\ddot{x}+c\dot{x}+kx=0$
$x(0)=x_0,\dot{x}(0)=v_0$

from data points of the form $(t_i,x_i)$ with $i=0,1,2,\dots ,n$.

Note. The Python program embedded below will generate simulated data for us.

Method of Parameter Estimation.

We will think of the parameters that we want to estimate, the $c$ and the $k$, as “living” in parameter space which consists of points $(c,k)$.

At each point $(c,k)$ in parameter space we can solve the differential equation

$\ddot{x}+c\dot{x}+kx=0$
$x(0)=x_0,\dot{x}(0)=v_0$

We write the solution as $x(c,k,t)$, instead of as just $x(t)$, to emphasize that the solution is a function of of the parameters $c$ and $k$.

We want to find the $(c,k)$ which minimizes the distance between $x(c,k,t)$ and the data. We measure this distance using the Sum of the Squared Errors (SSE).

We can think of the DE’s solution, $x(c,k,t)$, as predicting the position of the mass at time $t.$

We calculate $x(c,k,t)$ using the Euler method.

The data point $(t_i,x_i)$ is the actual or true position of the mass at time $t$. The $i^{th}$ error is the difference between where $x$ predicted the mass would be at time $t_i$, which is $x(c,k,t_i)$ and were the mass actually is at time $t_i$, which is $x_i$.

For each $i,i=0,1,2,\dots ,n$, the error is
$${ϵ}_i=x(c,k,t_i)–x_i$$
and the squared error is
$${ϵ}_i^2={(x(c,k,t_i)–x_i)}^2$$
Finally, to get the SSE, we sum up the squared errors
$$SSE(c,k)=\sum _{i=0}^n{(x(c,k,t_i)–x_i)}^2$$
Note. The SSE is a a function of $(c,k)$, hence the notation $SSE(c,k)$.

Keep in mind, each time the computer calculates the $SSE$ at $(c,k)$ the computer has run the Euler method algorithm to calculate the Euler points which approximate  $x(c,k,t)$. Then the computer uses interpolation to estimate each $x(c,k,t_i)$. For more about interpolation see the notes at the bottom of this page.

Note about SSE. For some general information on the SSE in the context of regression see:
https://mccarthymat150.commons.gc.cuny.edu/units-16/17-correlation-and-regression/

Next, we need an efficient algorithm to find the $(c,k)$ which minimizes $SSE(c,k)$. The algorithm we will use is called the Gradient Descent Algorithm (GDA).

Note about GDA. For some general information about the GDA see:
https://mccarthymat501.commons.gc.cuny.edu/gradient-descent-algorithm-gda/

The following Python 3 program  uses entirely numerical methods to estimate the parameters $(c,k).$ It solves the harmonic oscillator differential equation using Euler’s Method. It calculates the gradients used in the GDA numerically.

Outline of how the program works. You can accept the defaults, or you can chose values for $c$ and $k$ to be used in the differential equation model. Similarly, you can accept the defaults, or you can choose a different starting point $(c_0,k_0)$ for the GDA. Another important variable you can change is called

$\text{ GDAiterations }$

its default value is 200, meaning the GDA should iterate itself 200 times.

The program then uses the Euler method to numerically solve the differential equation for the chosen values of $c$ and $k$. It uses that solution to generate data with noise in it using a random number generator. The idea is to simulate the sorts of measurements you would get if you were working with a physical damped harmonic oscillator, like in a physics lab.

The program applies the GDA to the $SSE(c,k)$ to estimate the parameters $(c,k)$.

The program outputs the GDA estimate for $(c,k)$ (which it reports as c/m and k/m) as well as the SSE for that estimate and some additional information:

Then, if the option $\text{ Do3DandContourPlot = 1 }$ the program will output two graphs.

Note. The most time consuming part of the program is doing the calculations to produce these graphs. As these calculations are done, Python outputs what percent of the calculations are left to do.

The first is the SSE(c,k) surface with the GDA path superimposed on it.  The surface height and GDA path height are scaled by log 10. This is to prevent the huge values of the SSE at some locations from making it impossible to see, in the graph, the details at locations were the SEE is near its minimum.  See image below:

Note. All images were produced by the source code embedded at the bottom of this page.

Note that in Python, the above graph can be rotated interactively:

The second graph is the SSE(c,k) surface, again scaled by log 10, shown as a two dimensional contour map. The contours being the level curves of the surface. The GDA path is superimposed on it.  See image below:

Note that in Python, sections of the graph can be magnified using the magnifying glass tool. That tool is at the top of the graph window. See video and image below:

Video of Contour Map showing SSE(c,k) and the GDA path in (c, k) space.

In the above image, which is a zoomed in version of the previous image, we see the last points produced by the GDA. Notice how the GDA nicely gets us close to the minimum SSE (the center of that blue polygonal shape) but then starts to zigzag a little. That zigzagging  behavior is a common artifact of GDA. To counter that, I have our GDA output the $(c,k)$ estimate as the average of the last two GDA points.

After the optional 3D and contour graphs, the program will produce a graph showing the data points $(t_i,x_i)$ (red dots) together with the solution $(t,x(t))$ with x(t) being numerically solved by applying the Euler method to the initial choice of c, k, i.e., from the model (yellow line) and from the final GDA estimate of c, k (blue line).  See graph below.

The program will finally produce a video of the above graph, showing the evolution the solution x(c,k, t) as the GDA updates its estimates of (c, k). See video below.

Note. Unlike the static graph above, that the video does not graph $x(c,k,t)$ for the final estimate of  $(c,k)$. Recall that the final estimate of $(c,k)$ comes from the average of the final two GDA estimates for $(c,k)$. We take that average in order to compensate for the zigzagging that tends to happen at the end of the GDA.

In the video below, the up and down motion of the $x(c,k,t)$ (blue) starting at GDA iteration 170 is due to the GDA estimates of $(c,k)$ zigzagging. This zigzagging can be seen towards the end of the video of the contour map (embedded a few paragraphs above).

 

Some technical notes about the above Python implementation.

1. Linear Interpolation.  When the Euler method is used to numerically solve a differential equation what is actually produced is a finite set of discrete points $(t_j,x_j)$ called Euler points.  We convert the Euler points into a continuous function using Linear Interpolation. Basically, we connect the Euler points with straight line segments. See image below.
2. We actually use a modification of the classical GDA.
   What we use is an example of an adaptive GDA. We adapt the learning rate $\eta$ to the size of the gradient. If the length of the gradient vector $ng>1$ we replace $\eta$ with $\eta /ng$. This has the effect of insuring that the distance between GDA points is at most $\eta$. Note. $ng$ stands for “norm of the gradient” which is a fancy way to say the length of the gradient vector.
   Note. The GDA method is based on the fact that the negative of the gradient points in the direction of fast descent. However, this is a local property. We don’t know for how far we can go in that direction (in the direction of the negative of the gradient) and still be descending.  If we move too far in that direction, instead of descending, we might be ascending. This is like if you jump over a hole and instead of descending into it, you end up on higher ground. The entire purpose of the GDA is to go down into the hole (to reach a local minimum).