Notebook

HARMONIC OSCILLATORS¶

(A tutorial on modelling and simulation with Python)¶

1. Harmonic oscillator¶

An Harmonic Oscillator is a dynamical system that presents only a single stable attractor or equlibrium point for a particular time-dependent variable . The principal characteristic about harmonic oscillators is that they experience a restoring force with a magnitude proportional to the gap between the variable and its equilibrium point.

More formally, a harmonic oscilator can be presented as a Linear Ordinary Differential Equation (Linear ODE) stated as follows:

$a \frac{d^2y}{dt} = -cy$

Where $y = y(t)$ is the time-dependent variable, $a\frac{d^2y}{dt}$ is the restoring force experienced by $y$ , and $-cy$ is the magnitude of such restoring force proportional to the distance from equilibrium and in opposite direction. $a,c$ are constant coefficients

2. Solution of Hamonic oscillator¶

The name of Harmonic Oscillator comes because its solution consist of oscillating functions such as sine or cosine, more specifically:

$y(t) = A\textrm{sin} \left( \sqrt{\frac{c}{a}}t + \phi_1 \right) + B\textrm{cos} \left( \sqrt{\frac{c}{a}}t + \phi_2 \right)$

Where $\phi_1$ and $\phi_2$ are phase-shift constants and $A,B$ are amplitude constants.

Demonstration:¶

$y(t) = A\textrm{sin} \left( \sqrt{\frac{c}{a}}t + \phi_1 \right) + B\textrm{cos} \left( \sqrt{\frac{c}{a}}t + \phi_2 \right)$

$\frac{dy(t)}{dt} = \sqrt{\frac{c}{a}}A\textrm{cos} \left( \sqrt{\frac{c}{a}}t + \phi_1 \right) - \sqrt{\frac{c}{a}}B\textrm{sin} \left( \sqrt{\frac{c}{a}}t + \phi_2 \right)$

$\frac{d^2y(t)}{dt^2} = -\left(\sqrt{\frac{c}{a}}\right)^2A\textrm{sin} \left( \sqrt{\frac{c}{a}}t + \phi_1 \right) - \left(\sqrt{\frac{c}{a}}\right)^2B\textrm{cos} \left( \sqrt{\frac{c}{a}}t + \phi_2 \right)$

$\frac{d^2y(t)}{dt^2} = -\frac{c}{a}\left( A\textrm{sin} \left( \sqrt{\frac{c}{a}}t + \phi_1 \right) + B\textrm{cos} \left( \sqrt{\frac{c}{a}}t + \phi_2 \right) \right)$

$\frac{d^2y(t)}{dt^2} = -\frac{c}{a} y(t)$

$a\frac{d^2y(t)}{dt^2} = -c y(t)$

Python Example:¶

In [1]:

import numpy as np
import matplotlib.pyplot as plt

c = 1
a = 1
A = 1
B = 1
phi_1 = 0 
phi_2 = 0

t = np.linspace(0,2*np.pi,1000); #Timer Vector
y = A*np.sin(np.sqrt(c/a)*t+phi_1) + B*np.cos(np.sqrt(c/a) + phi_2) #y(t)
d2y_dt2 = np.diff(np.diff(y))/(0.005)**2 #second derivative of y(t)

#Plot results
plt.grid(True)
plt.plot(t,y,label='$y(t)$')
plt.plot(t[2:],d2y_dt2,label='$\ddot{y}(t)$')
plt.legend()
plt.show()

3. The damped Harmonic Oscillator¶

The Harmonic Oscillator looks nice for waving systems, such as pendulums and mass-spring, where the system exchanges potential and kinetic energy harmoniously in an endless waving movement.

However, the reality does not behavior like that. In the real world, harmonic systems have energy losses due to the second-law of thermodynamic that states that it is impossible to convert from one form of energy to another with 100% of efficiency. Additionally, the faster is the process of conversion of energy, the more energy losses.

Such fastness can be interpreted as the rate-of-change of the a variable in time, i.e. its first-derivative. So the Harmonic system equation:

$a \frac{d^2y}{dt^2} = -cy$

becomes

$a \frac{d^2y}{dt^2} = -cy - b \frac{dy}{dt}$

$a \frac{d^2y}{dt^2} + b\frac{dy}{dt} + cy = 0$

Where $b$ is known as damping factor. The higher the value of $b$ , the more the energy losses over time.

4. Numerical methods for the solution of the damped Harmonic system¶

Despite the solution of damped Harmonic system can be found by hand, we want to introduce numerical methods to solve dynamical equations using Python. We will use this knowledge later to solve systems that can not be solved analytically.

Linear ODE's can be expressed as finite-difference equations to be solved numerically (using computational algorithms). For this case, follow these steps:

Use an auxiliar variable $z$ to express the original second-order equation as a first-order equation system

$a \frac{d^2y}{dt^2} + b\frac{dy}{dt} + cy = 0,~~~~~~ z = \frac{dy}{dt}$

Replacing $z$ in the first equation we obtain the equation system

$a \frac{dz}{dt} + bz + cy = 0 ~~~~~\textrm{(1)}$

$\frac{dy}{dt}-z = 0 ~~~~~~~~~~~~~~~~~\textrm{(2)}$

$~$
Express all derivatives as finite-differences using Euler's method.

$a \left(\frac{z(t) - z(t-\Delta t)}{\Delta t }\right) + bz + cy = 0 ~~~~~\textrm{(1)}$

$\frac{y(t)-y(t-\Delta t)}{\Delta t}-z = 0 ~~~~~~~~~~~~~~~~~~~~~~\textrm{(2)}$

$~$
Bring the system to the discrete-time domain

$a \left(\frac{z[n] - z[n-1]}{\Delta t }\right) + bz[n] + cy[n] = 0 ~~~~~\textrm{(1)}$

$\frac{y[n]-y[n-1]}{\Delta t}-z[n] = 0 ~~~~~~~~~~~~~~~~~~~~~~~~~~~\textrm{(2)}$

$~$
Clear $y[n]$ in (2)

$a \left(\frac{z[n] - z[n-1]}{\Delta t }\right) + bz[n] + cy[n] = 0 ~~~\textrm{(1)}$

$y[n] = y[n-1] + \Delta t z[n]~~~~~~~~~~~~~~~~~~~~~~~~~~~~\textrm{(2)}$

$~$
Replace (2) in (1)

$a \left(\frac{z[n] - z[n-1]}{\Delta t }\right) + bz[n] + c \left(y[n-1] + \Delta t z[n] \right) = 0 ~~~\textrm{(1)}$

$y[n] = y[n-1] + \Delta t z[n]~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~\textrm{(2)}$

$~$
Clear z[n] in (1)

$z[n] = \left( \frac{a}{\Delta t}z[n-1] - cy[n-1] \right) / \left( \frac{a}{\Delta t} + b +c\Delta t \right) ~~~~~~~(1)$

$y[n] = y[n-1] + \Delta t z[n]~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~\textrm{(2)}$

$~$

Now, we can find the an aproximation of the solution $y[n]$ and its first derivative $z[n]$ using a logistic map algorithm that depends of certain initial conditions $y[0]$ and $z[0]$ and a sampling time $\Delta t$ .

Implememtation in Python¶

In [2]:

# Coefficients setup
points = 5000
a = 100.
b = 10.
c = 100.

# Time vector for the solution y[n] filled with zeros.
y = np.zeros((points,1))

# Initial conditions and sampling time
dt = 0.01;
y[0] = 50; #y[0]
dydt = 0;  #z[0]

#Simulation
for n in xrange(1,points):
	dydt = ((a/dt)*(dydt) - c*y[n-1]) / ((a/dt) + b + c*dt)
	y[n] = y[n-1] + dt*dydt

#Plot results
plt.grid(True)    
plt.plot(np.linspace(0,dt*points,points),y,label='y(t)')
plt.legend()
plt.show()

5. Response of Harmonic systems to external forces¶

Until now, we have analyzed the behavior of Harmonic systems as a closed system evolving in the time only due to the action of its intrinsic variables and a couple of initial conditions.

However, in the reality, harmonic systems may be affected or influenced by external forces, so, the response of the system changes significantly. (More if the external forces have a considerable magnitude).

To involve external forces in the harmonic system, just define the original equation of damped Harmonic oscillator as follows:

$a \frac{d^2y}{dt^2} + b\frac{dy}{dt} + cy = F(t)$

Where $F(t)$ is the external force changing in time that may be a constant function, or maybe periodic, linear, non-linear, even non-analytic.

Now is when numerical methods become a better choice over analytical methods, because when $F(t)$ is non-linear, the solution of harmonic system is pretty hard to obtain. Even worse, if $F(t)$ can not be expressed as linear combination of analytical functions (like the noise), the solution is impossible to get by hand.

In contrast, numerical methods does not care about the shape of $F(t)$ , if it can be discretized, the algorithm just evolves in time steps using the information of $F[n]$ and the previous values of $y$ and $z$ , i.e. $y[n-1]$ and $z[n-1]$ .

$a \frac{d^2y}{dt^2} + b\frac{dy}{dt} + cy = F(t)$

becomes

$z[n] = \left(F[n] + \frac{a}{\Delta t}z[n-1] - cy[n-1] \right) / \left( \frac{a}{\Delta t} + b +c\Delta t \right) ~~~~~~~(1)$

$y[n] = y[n-1] + \Delta t z[n]~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~\textrm{(2)}$

$~$

Impulse response and step response of harmonic system in Python¶

In [3]:

# Coefficients setup
points = 5000
a = 100.
b = 10.
c = 100.

# Time vector for the solution y[n] filled with zeros.
y = np.zeros((points,1))

# Initial conditions and sampling time
dt = 0.01;
y[0] = 50; #y[0]
dydt = 0;  #z[0]

# Impulses
f = np.zeros((points,1))
for i in xrange(2):
    f[i*2500] = 50
    
#Simulation
for n in xrange(1,points):
	dydt = f[n] + ((a/dt)*(dydt) - c*y[n-1]) / ((a/dt) + b + c*dt)
	y[n] = y[n-1] + dt*dydt

#Plot results
plt.grid(True)    
plt.plot(np.linspace(0,dt*points,points),y,label='y(t)')
plt.plot(np.linspace(0,dt*points,points),f,label='F(t)')
plt.legend()
plt.show()

In [4]:

# Coefficients setup
points = 5000
a = 100.
b = 10.
c = 100.

# Time vector for the solution y[n] filled with zeros.
y = np.zeros((points,1))

# Initial conditions and sampling time
dt = 0.01;
y[0] = 50; #y[0]
dydt = 0;  #z[0]

# Step
f = np.zeros((points,1))
for i in xrange(2):
    if (i+1)%2  == 0:
        f[i*2500:(i+1)*2500] = 0.5;

#Simulation
for n in xrange(1,points):
	dydt = f[n] + ((a/dt)*(dydt) - c*y[n-1]) / ((a/dt) + b + c*dt)
	y[n] = y[n-1] + dt*dydt

#Plot results
plt.grid(True)    
plt.plot(np.linspace(0,dt*points,points),y,label='y(t)')
plt.plot(np.linspace(0,dt*points,points),f*100,label='F(t)')
plt.legend()
plt.show()

6. Applications: RLC circuit as a band-pass filter¶

For a series RLC circuit, the voltage of the source $V_s(t)$ is divided between the resistor $R$ , the capacitor $C$ and the inductor $L$ following the Kirchhoff's voltage law:

$V_L(t) + V_R(t) + V_C(t) = V_s(t)$

$~$

$L\frac{dI(t)}{dt} + RI(t) + \frac{Q(t)}{C} = V_s(t), ~~~~ \textrm{where} ~~~~ I(t) = \frac{dQ(t)}{dt}$

Where, $I(t)$ is the common current in amperes for the series loop, and $Q(t)$ is the electric charge in coulombs. The rate-of-change of the charge along the time is equal to the current.

This formulation looks pretty similar to that described in section 5. In fact, $a$ becomes $L$ , the damping coefficient $b$ is the resistance $R$ , $1/C$ is equal to $c$ , the dynamical variable is $Q(t)$ and the fastness of the variable is defined by $I(t)$ .

Problem¶

We want to see the response-voltage on the resistor for a periodic signal as $Vs(t)$ with fundamental frequency $f_0=500$ Hz

Model the equation system $L\frac{dI(t)}{dt} + RI(t) + \frac{Q(t)}{C} = V_s(t) ~~ \textrm{(1)}$ $I(t) = \frac{dQ(t)}{dt} ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~\textrm{(2)}$ $~$
Simulate the model to find the response of the system in terms of $Q(t)$ , $I(t)$ , and $V_R(t) = RI(t)$

In [5]:

# Coefficients setup
points = 10000 #0.1 seconds with dt = 0.00001
L = 10e-3 #10 mH
R = 30   #30  ohm
C = 10e-6 #10 uF

# Time vector for the solution Q[n] and I[n] filled with zeros.
Q = np.zeros((points,1))
I = np.zeros((points,1))


# Initial conditions and sampling time
dt = 0.00001; # 10 us
Q[0] = 0;  #y[0]
I[0] = 0;  #z[0]

# Source voltage
t = linspace(0,0.1,10000)
f = 500 # natural frequency 500Hz 
Vs = 10*np.cos(2*np.pi*f*t)


#Simulation
for n in xrange(1,points):
	I[n] = (Vs[n] + (L/dt)*(I[n-1]) - (1/C)*Q[n-1]) / ((L/dt) + R + (1/C)*dt)
	Q[n] = Q[n-1] + dt*I[n]

#Voltage on the resistor
VR = R*I

#Plot results
plt.figure()
plt.grid(True)    
plt.ylim(-15,15)
plt.plot(np.linspace(0,dt*points,points),Vs,'b-',label='$V_s(t)$')
plt.legend()
plt.show()

plt.grid(True)    
plt.ylim(-15,15)
plt.plot(np.linspace(0,dt*points,points),VR,'g-',label='$V_R(t)$')
plt.legend()
plt.show()

This RLC circuit example has a natural frequency $f_0 = 500$ Hz. That is, the voltage on the resistor will be maximum at 500Hz. By modifying the values of $L$ and $C$ the value of $f_0$ will be different. Additionally, by adjusting R, we affect the bandwidth of the circuit. A circuit with high $R$ becomes very selective, supressing frequencies out of $f_0$ . In contrast, low R values just supress frequencies that are far away from $f_0$ .

*That is how radio-frequencies (RF) communications work!*, a carrier signal with fundamental frequency $f_0$ is sent to the air and the receptor will be able to receive the signal only if is adjusted to receive that fundamental frequency.

You are invited to modify $L$ , $C$ , $R$ and $f$ in the code to see different responses of the system

Now is your turn!¶

7. Problem: Modelling of the response for a vehicle's suspension¶

A shock-absorber is nothing but a mass-spring system that can be modeled as an Harmonic Oscillator. However, in this case we do not want to see the natural response of the shock-absorber as a closed system but its response to an applied external force, e.g. the bumps on the road.

A good shock-absorber must be rigid if the force has no strong variations over time, also if the force is applied in a very short time. On the other hand, the shock-absorber must compress itself if the force varies considerably in the order of tenths of second.

Example

In [62]:

def gaussian(x, mu, sig):
    return (1/np.sqrt(2*np.pi*sig**2))*np.exp(-np.power(x - mu, 2.) / (2 * np.power(sig, 2.)))

t = linspace(0,4,1000);

plt.grid(True)
plt.ylim(0,150);
width = 0.04;
plt.plot(t,10*gaussian(t,2,width),'g-',linewidth=2);

width = 0.1;
plt.plot(t,10*gaussian(t,2,width),'r-',linewidth=2,label='must be absorbed');

width = 1;
plt.plot(t,10*gaussian(t,2,width),'g-',linewidth=2,label='must not be absorbed');

plt.legend()

Out[62]:

<matplotlib.legend.Legend at 0x7f7ab4af5110>

Problem For the applied forces described above, we want to adjust the values of the spring constant $k$ and the damping factor $c$ for a mass $m = 30$ kg for a vehicle suspension such that the system oscillates for medium-duration forces and become rigid otherwise.

Tips

Model the problem as an harmonic oscillator
Use a similar code to those used in this tutorial to simulate the system
Adjust parameters to find the adequate response

You reached the end of this tutorial¶