In [1]:

```
%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
from scipy.linalg import solve
```

Beginning in 1595, Galileo began a series of experiments on projectile motion. In one such experiment, he let a ball roll down an inclined ramp right off of a table and measured it's vertical distance travelled. He did this for several ramp heights and thus obtained a relationship between the height at which he released the ball and the vertical distance that the ball traveled. I've got the actual data stored on my website! We can load and view it directly with Python like so:

In [2]:

```
df = pd.read_csv('https://marksmath.org/data/galileos_ramp.csv')
df
```

Out[2]:

The unit of distance is the "punti", according to the source:

We can also plot the data:

In [3]:

```
plt.plot(df.D,df.H, '.')
```

Out[3]:

While it rather flips the role of what we might call *explanatory* and *response*, it does look like height released is related to distance traveled via a quadratic - that is,
$$H = a D^2 + b D + c.$$
In order to fit the data, we could solve the least squares problem
$$
\left(\begin{matrix}
d_1^2 & d_1 & 1 \\
d_2^2 & d_2 & 1 \\
\vdots & \vdots & \vdots \\
d_n^2 & d_n & 1
\end{matrix}\right)
\left(\begin{matrix}a \\ b \\ c \end{matrix} \right) \approx
\left(\begin{matrix}
h_1 \\ h_2 \\ \vdots \\ h_n
\end{matrix}\right).
$$

Of course, using the normal equations $A^T A \vec{x} = A^T \vec{b}$, this boils down to solving the linear system

$$ \left(\begin{matrix} d_1^2 & d_2^2 & \cdots & d_n^2 \\ d_1 & d_2 & \cdots & d_n \\ 1 & 1 & \cdots & 1 \end{matrix}\right) \left(\begin{matrix} d_1^2 & d_1 & 1 \\ d_2^2 & d_2 & 1 \\ \vdots & \vdots & \vdots \\ d_n^2 & d_n & 1 \end{matrix}\right) \left(\begin{matrix}a \\ b \\ c \end{matrix} \right) = \left(\begin{matrix} d_1^2 & d_2^2 & \cdots & d_n^2 \\ d_1 & d_2 & \cdots & d_n \\ 1 & 1 & \cdots & 1 \end{matrix}\right) \left(\begin{matrix} h_1 \\ h_2 \\ \vdots \\ h_n \end{matrix}\right). $$Here's how to set up the matrix containing the values of $D$ - the matrix traditionally called $A$ in the normal equations:

In [4]:

```
A = np.array([df.D**2, df.D, df.D**0]).transpose()
A
```

Out[4]:

We now solve the normal equations:

In [5]:

```
a,b,c = solve(A.transpose().dot(A), A.transpose().dot(df.H))
a,b,c
```

Out[5]:

Thus we expect $$H \approx 0.00818073977 D^2 - 4.00484297 D + 600.05.$$ Let's define and plot that function together with data:

In [6]:

```
def f(x): return a*x**2 + b*x + c
xs = np.linspace(250,600)
ys = f(xs)
plt.plot(xs,ys)
plt.plot(df.D,df.H, '.')
```

Out[6]:

Looks pretty good!