Numerical integration of coupled first-order ODEs¶
Uses odeint() function from the integrate submodule of scipy.
Marty Ligare, August 2020
Damped harmonic oscillator¶
The equation of motion for a damped harmonic oscillator is $$ \frac{d^2x}{dt^2} = -\omega_0^2\, x - \gamma \frac{dx}{dt}. $$ Working with dimensionless variables in which time is measured in units of $\omega_0^{-1}$ this equation of motion can be rewritten as $$ \frac{d^2x}{dt^2} = -x - b\frac{dx}{dt}, $$ where $b = \gamma/\omega_0$.
This second order differential equation can be written as two coupled first-order equations:
\begin{eqnarray*} \frac{dx}{dt} &=& v\\ \frac{dv}{dt} &=& -x - bv \end{eqnarray*}The function below returns the RHS of these differential equations; the lines above the return relates the elements of the array u[i] with variables that have more physical meaning.
import numpy as np
from scipy import integrate
import matplotlib as mpl
import matplotlib.pyplot as plt
# Following is an Ipython magic command that puts figures in the notebook.
%matplotlib notebook
# Following sets up LateX fonts
#mpl.rc('font', **{'family': 'serif', 'serif': ['Computer Modern']})
#mpl.rc('text', usetex=True)
# M.L. modification of matplotlib defaults
# Changes can also be put in matplotlibrc file,
# or effected using mpl.rcParams[]
plt.style.use('classic')
plt.rc('figure', figsize = (6, 4.5)) # Reduces overall size of figures
plt.rc('axes', labelsize=14, titlesize=14)
plt.rc('xtick', labelsize=12)
plt.rc('ytick', labelsize=12)
plt.rc('figure', autolayout = True) # Adjusts supblot parameters for new size
def damped_osc(u,t):
'''Function returning first-order derivatives of the dependent
quantities x and v. The values of x and v are passed in the
two-element array u, where u[0] and u[1] correspond to x and v
respectively.'''
x = u[0]
v = u[1]
return (v,-x-b*v)
Initial conditions, and parameter(s)¶
x0 = 2 # initial position
v0 = 0 # initial velocity
b = 0.4 # camping parameter
u0 = np.array([x0,v0])
In following cell:
- Choose points for output.
- Integrate.
Note:odeintreturns an array:
[[x_0 v_0],
[x_1 v_1],
[x_2 v_2], ...]
To get single list forxand single list forvwe need the transpose of the returned array.
(Could also keep return as a single array if that's more useful down the road.)
t = np.linspace(0,20,201) # NOTE: The points selected for plotting are
# not the points used for the numerical
# evalution.
x, v = integrate.odeint(damped_osc,u0,t).T
plt.figure()
plt.plot(t,x)
plt.axhline(0)
plt.xlabel("$t$")
plt.ylabel("$x(t)$");
plt.figure()
plt.plot(t,v);
plt.axhline(0)
plt.xlabel('$t$')
plt.ylabel('$v(t)$');
plt.figure()
plt.xlabel("$x$")
plt.ylabel("$v$")
plt.title("phase space",fontsize=14)
plt.plot(x,v);
Version Information¶
version_information is from J.R. Johansson (jrjohansson at gmail.com); see Introduction to scientific computing with Python for more information and instructions for package installation.
version_information is installed on the linux network at Bucknell
%load_ext version_information
version_information numpy, scipy, matplotlib