import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit

plt.rcParams["figure.figsize"] = (10,6)
plt.rcParams.update({'font.size': 14})

1Basilar Membrane¶

We learn that the basilar membrane is a membrane that gains mass and loses stiffness as it extends further from the oval window toward the apex. We model it as a spring-mass system with varying mass and spring stiffness. Let’s look at how those variables affect vibrational frequency.

According to Hooke’s Law

F = -kx

(1)

saying a spring with stiffness $k$ will pull back against an applied force. The force required to extend or compress a spring a distance $x$ is $-kx$ . We can apply a force using gravity to verify that a spring is “Hookian”. In this case, adding a mass to a spring and allowing gravity to pull the mass downward will create an upward force in the spring to balance the gravitational force.

\sum F = kx-mg = 0

(2)

mg = kx

(3)

I’ve made the spring force $kx$ positive for upward and the weight $mg$ downward or negative. If we plot $mg$ vs. $x$ , the slope will be $k$ . Let’s add mass and measure the extension of the spring. Enter your quantities below.

g = 9.8 #m/s^2
m = np.array([0.05,0.06,0.07,0.08])#array of masses
mg = m*g
x = np.array([0.05, 0.06, 0.07, 0.08])#array of extension lengths in meters
x_unc = np.array([0.005,0.005,0.005, 0.005])#estimate your uncertainty in measuring x
mg_unc = x_unc/x*mg#propagate the relative uncertainty in x to mg


def f_line(x, m, b):
    return m*x+b

fit_params, fit_cov =curve_fit(f_line, x, mg, p0=(1,1), sigma=mg_unc, absolute_sigma=True)
print(fit_params)

y_fit = f_line(x, fit_params[0], fit_params[1])

plt.errorbar(x, mg, yerr=mg_unc, fmt='ok', ecolor=None, elinewidth=1.0, capsize=1.0, capthick=1.0)
plt.plot(x, y_fit, '-k', label='fit')
plt.ylabel('mg (N)')
plt.xlabel('x (m)')
plt.show()

1.1Solving the Equation of Motion¶

If we attach a mass to a spring or simply have something massive that is elastic in a “Hookian” way, we could stretch it by pulling on it. The force we apply to stretch it satifies the equation

F_{pull} = -kx

(4)

This is a $F_{net}=0$ situation like the one described above. If we let go, the elastic material will have a net force that comes from Hooke’s Law.

F_{net} = -kx

(5)

Since there is a net force, there is acceleration.

F_{net} = -kx = ma

(6)

and since the acceleration is along $x$ , we can rewrite this as an “equation of motion”

-kx = \frac{d^2x}{dt^2}

(7)

This is a commonly recurring differential equation in physics, and we know the solutions are

x(t) = \begin{cases} A\cos\omega t\\ A\sin\omega t \end{cases}

(8)

where ω is the frequency the spring-mass will oscillate. We can solve the equation of motion with each of these and see that

\omega = \sqrt{\frac{k}{m}}

(9)

and this can be rewritten for the period of oscillation $T= 1/f = 2\pi/\omega$ .

T = 2\pi\sqrt{\frac{m}{k}}

(10)

1.2Frequency vs. mass¶

Now, let’s look at how the period changes with mass by doing an experiment. We’ll record the period of oscillation for various masses. Use a stopwatch to measure the period (top-to-top cycle for example). It may be more accurate to measure 3 to 5 periods and divide by the number of periods measured. Enter your results below, and fit the results with a power law function. To do the fit, enter the left hand side of the function

T = ax^b

(11)

From the equation for period in terms of mass and spring constant, what are you expecting for $a$ and $b$ ?

m = np.array([0.05, 0.06, 0.07, 0.08, 0.09, 0.10])#array of masses
T = np.array([0.05, 0.06, 0.07, 0.08,  0.09, 0.10])#array of periods for different masses
T_unc = np.array([0.005, 0.005, 0.005, 0.005, 0.005, 0.005])#estimate your uncertainty in measuring T

def f_pow(x, a, b):
    return 

fit_params, fit_cov =curve_fit(f_pow, m, T, p0=(1,-1), sigma=T_unc, absolute_sigma=True)
print(fit_params)

y_fit = f_pow(m, fit_params[0], fit_params[1])

plt.errorbar(m, w, yerr=w_unc, fmt='ok', ecolor=None, elinewidth=1.0, capsize=1.0, capthick=1.0)
plt.plot(m, y_fit, '-k', label='fit')
plt.ylabel('T (s)')
plt.xlabel('mass (kg)')
plt.show()