This example shows how the Simple Harmonic Oscillator model can be used.
A model for a particle undergoing Newtonian dynamics that experiences a force in proportion to its displacement from an equilibrium position, and, in addition, a friction force. The motion of the particle can be determined by solving a second order ordinary differential equation (from Newton's $F = ma$):
$$\frac{d^2y}{dt^2} = -y(t) - \theta \frac{dy(t)}{dt}.$$where $y(t)$ is the particle's displacement and $\theta$ is the frictional force.
In [1]:
import pints
import pints.toy
import matplotlib.pyplot as plt
import numpy as np
model = pints.toy.SimpleHarmonicOscillatorModel()
Parameters are given in the order $(y(0), dy/dt(0), \theta)$. Here, we see that since $\theta > 0$, that the oscillations of the particle decay exponentially over time.
In [7]:
times = np.linspace(0, 50, 1000)
parameters = model.suggested_parameters()
values = model.simulate(parameters, times)
plt.figure(figsize=(15,2))
plt.xlabel('t')
plt.ylabel(r'$y$ (Displacement)')
plt.plot(times, values)
plt.show()
Substituting an exponential solution of the form: $y(t) = Ae^{\lambda t}$ into the governing ODE, we obtain: $\lambda^2 + \theta \lambda + 1=0$, which has solutions:
$$\lambda = (-\theta \pm \sqrt{\theta^2 - 4})/2.$$As we can see above, if $\theta^2 < 4$, i.e. if $-2<\theta<2$, $\lambda$ has an imaginary part, which causes the solution to oscillate sinusoidally whilst decaying to zero.
If $\theta = 2$, we get critically dampled dynamics where the displacement decays exponentially to zero, rather than oscillatory motion.
In [8]:
parameters = model.suggested_parameters()
values = model.simulate([1, 0, 2], times)
plt.figure(figsize=(15,2))
plt.xlabel('t')
plt.ylabel(r'$y$ (Displacement)')
plt.plot(times, values)
plt.show()
If $\theta > 2$, we get overdamped dynamics, with an increased rate of decay to zero.
In [10]:
parameters = model.suggested_parameters()
values = model.simulate([1, 0, 5], times)
plt.figure(figsize=(15,2))
plt.xlabel('t')
plt.ylabel(r'$y$ (Displacement)')
plt.plot(times, values)
plt.show()
This model also provides sensitivities: derivatives $\frac{\partial y}{\partial p}$ of the output $y$ with respect to the parameters $p$.
In [5]:
values, sensitivities = model.simulateS1(parameters, times)
We can plot these sensitivities, to see where the model is sensitive to each of the parameters:
In [12]:
plt.figure(figsize=(15,7))
plt.subplot(3, 1, 1)
plt.ylabel(r'$\partial y/\partial y(0)$')
plt.plot(times, sensitivities[:, 0])
plt.subplot(3, 1, 2)
plt.xlabel('t')
plt.ylabel(r'$\partial y/\partial \dot y(0)$')
plt.plot(times, sensitivities[:, 1])
plt.subplot(3, 1, 3)
plt.xlabel('t')
plt.ylabel(r'$\partial y/\partial \theta$')
plt.plot(times, sensitivities[:, 1])
plt.show()