Numerical integration of coupled first-order ODEs¶
Damped harmonic oscillator¶
The equation of motion for a damped harmonic oscillator is ta $$ \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 scipy as sp
from scipy import integrate # not included in basic scipy
import matplotlib as mpl # As of July 2017 Bucknell computers use v. 2.x
import matplotlib.pyplot as plt
# Following is an Ipython magic command that puts figures in the notebook.
# For figures in separate windows, comment out following line and uncomment
# the next line
# Must come before defaults are changed.
%matplotlib notebook
#%matplotlib
# As of Aug. 2017 reverting to 1.x defaults.
# In 2.x text.ustex requires dvipng, texlive-latex-extra, and texlive-fonts-recommended,
# which don't seem to be universal
# See https://stackoverflow.com/questions/38906356/error-running-matplotlib-in-latex-type1cm?
mpl.style.use('classic')
# M.L. modifications of matplotlib defaults using syntax of v.2.0
# More info at http://matplotlib.org/2.0.0/users/deflt_style_changes.html
# Changes can also be put in matplotlibrc file, or effected using mpl.rcParams[]
plt.rc('figure', figsize = (6, 4.5)) # Reduces overall size of figures
plt.rc('axes', labelsize=16, titlesize=14)
plt.rc('figure', autolayout = True) # Adjusts supblot parameters for new size
# Function returning derivatives of the dependent quantities u[0]
# and u[1], or more physically, x and v.
def damped_osc(u,t,b):
x = u[0]
v = u[1]
return (v,-x-b*v)
Initial conditions, and parameter(s)¶
x0 = 2
v0 = 0
b = 0.4 # Damping parameter
u0 = sp.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 = sp.linspace(0,20,201) # NOTE: The points selected for plotting are
# not the points used for the numerical
# evalution.
x, v = sp.integrate.odeint(damped_osc,u0,t,args=(b,)).T
plt.figure(1)
plt.plot(t,x);
plt.plot(t,v);
plt.axhline(0)
plt.xlabel("$t$")
plt.ylabel("$x(t)$ and $v(t)$")
plt.figure(2)
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:
http://nbviewer.jupyter.org/github/jrjohansson/scientific-python-lectures/blob/master/Lecture-0-Scientific-Computing-with-Python.ipynb
for more information and instructions for package installation.
If version_information has been installed system wide (as it has been on linuxremotes), continue with next cell as written. If not, comment out top line in next cell and uncomment the second line.
%load_ext version_information
#%install_ext http://raw.github.com/jrjohansson/version_information/master/version_information.py
version_information scipy, matplotlib