In this notebook, we'll be working on a classic problem: solving the harmonic oscillator equation.
where $\omega_0^2 = \frac{k}{m}$. The above equation is the harmonic oscillator model equation with a source term.
Since we have a second derivative in time of the position $x$, we need to specify two things: the initial position and the initial speed (i.e. time derivative) of the oscillator. These are denoted as:
$$ \begin{align} x(t = 0) = x_0 \\ \dot{x}(t=0) = \dot{x}_0 \end{align} $$In this notebook, we will explore three options for solving the evolution problem of this harmonic oscillator:
sympy
scipy
builtin toolsTo solve this equation analytically, we will use sympy. Sympy provides an ordinary differential equation (ODE) module for these problems: http://docs.sympy.org/dev/modules/solvers/ode.html.
First, we define our symbols and function:
import sympy
sympy.init_printing()
m, k, x_0, xdot_0, omega_0, t, omega = sympy.symbols('m, k, x_0, xdot_0, omega_0, t, omega')
x = sympy.Function('x')
We can then use dsolve
, which deals with differential equations:
sol = sympy.dsolve(sympy.Derivative(x(t), t, 2) + omega_0**2 * x(t) - sympy.sin(omega * t))
sol
And define our initial conditions and solve for them:
ics = [sympy.Eq(sol.args[1].subs(t, 0), x_0), sympy.Eq(sol.args[1].diff(t).subs(t, 0), xdot_0)]
ics
solved_ics = sympy.solve(ics)
solved_ics
full_sol = sol.subs(solved_ics[0])
full_sol
sympy.simplify(full_sol)
full_sol.subs({x_0:0, xdot_0:0})
The above equation was obtained by sympy and contains the solution to our problem. We can notice that it features complex exponentials, hinting at the oscillatory functions we already expect from our physical knowledge of harmonic oscillator.
Let's plot the solution for two different initial condition sets:
We'll use a value of $\omega_0$ equal to 2 (it's often helpful to not choose a value of 1 for constants to spot mistakes early).
case1 = sympy.simplify(full_sol.subs({x_0:1, xdot_0:0, omega_0:2, omega:1}))
case1
sympy.plot(case1.rhs, ylim=(-1.5, 1.5), xlim=(0, 10))
<sympy.plotting.plot.Plot at 0x24e9b79b4c8>
Let's look at our second solution.
case2 = sympy.simplify(full_sol.subs({x_0:0, xdot_0:1, omega_0:2, omega:1}))
case2
sympy.plot((case2.rhs, (t, 0, 20)), xlim=(0, 20), ylim=(-1, 1))
<sympy.plotting.plot.Plot at 0x24e9d3d1a48>
What about the limit case when we hit omega_0 ?
case2 = sympy.simplify(full_sol.subs({x_0:0, xdot_0:1, omega_0:2, omega:2}))
case2
cases = {}
for omega_value in [1, 1.4, 1.8, 1.9, 1.99]:
cases[omega_value] = make_solution(omega_value)
import numpy as np
import matplotlib.pyplot as plt
plt.style.use('bmh')
def make_plot(tmax=20):
t_values = np.linspace(0, tmax, num=1000)
fig, ax = plt.subplots(figsize=(15, 5))
for omega_value, case in cases.items():
func = sympy.lambdify(t, case.rhs, 'numpy')
ax.plot(t_values, func(t_values).real, label=r'$\omega$: {:.2f}'.format(omega_value))
ax.legend()
make_plot(tmax=20)
make_plot(tmax=100)
make_plot(tmax=1000)
cases = {}
for omega_value in [4 - 1, 4 - 1.4, 4 - 1.8, 4 - 1.9, 4 - 1.99]:
cases[omega_value] = make_solution(omega_value)
make_plot(tmax=20)
make_plot(tmax=100)
make_plot(tmax=1000)