{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Numerical integration of coupled first-order ODEs"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"#### Damped harmonic oscillator\n",
"The equation of motion for a damped harmonic oscillator is ta\n",
"$$ \\frac{d^2x}{dt^2} = -\\omega_0^2 x - \\gamma \\frac{dx}{dt}. $$\n",
"Working with dimensionless variables in which time is measured in \n",
"units of $\\omega_0^{-1}$ this equation of motion can be rewritten as \n",
"$$ \\frac{d^2x}{dt^2} = -x - b\\frac{dx}{dt}, $$\n",
"where $b = \\gamma/\\omega_0$.\n",
"\n",
"This second order differential equation can be written as two coupled first-order equations:\n",
"\n",
"\\begin{eqnarray*}\n",
"\\frac{dx}{dt} &=& v\\\\\n",
"\\frac{dv}{dt} &=& -x - bv\n",
"\\end{eqnarray*}\n",
"\n",
"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."
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [],
"source": [
"import scipy as sp\n",
"from scipy import integrate # not included in basic scipy\n",
"\n",
"import matplotlib as mpl # As of July 2017 Bucknell computers use v. 2.x \n",
"import matplotlib.pyplot as plt\n",
"\n",
"# Following is an Ipython magic command that puts figures in the notebook.\n",
"# For figures in separate windows, comment out following line and uncomment\n",
"# the next line\n",
"# Must come before defaults are changed.\n",
"%matplotlib notebook\n",
"#%matplotlib\n",
"\n",
"# As of Aug. 2017 reverting to 1.x defaults.\n",
"# In 2.x text.ustex requires dvipng, texlive-latex-extra, and texlive-fonts-recommended, \n",
"# which don't seem to be universal\n",
"# See https://stackoverflow.com/questions/38906356/error-running-matplotlib-in-latex-type1cm?\n",
"mpl.style.use('classic')\n",
" \n",
"# M.L. modifications of matplotlib defaults using syntax of v.2.0 \n",
"# More info at http://matplotlib.org/2.0.0/users/deflt_style_changes.html\n",
"# Changes can also be put in matplotlibrc file, or effected using mpl.rcParams[]\n",
"plt.rc('figure', figsize = (6, 4.5)) # Reduces overall size of figures\n",
"plt.rc('axes', labelsize=16, titlesize=14)\n",
"plt.rc('figure', autolayout = True) # Adjusts supblot parameters for new size"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"# Function returning derivatives of the dependent quantities u[0] \n",
"# and u[1], or more physically, x and v.\n",
"def damped_osc(u,t,b):\n",
" x = u[0]\n",
" v = u[1]\n",
" return (v,-x-b*v)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"#### Initial conditions, and parameter(s)"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"x0 = 2\n",
"v0 = 0\n",
"b = 0.4 # Damping parameter\n",
"u0 = sp.array([x0,v0])"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"In following cell:\n",
"- Choose points for output.
\n",
"- Integrate.
\n",
"Note: `odeint` returns an array:
\n",
" `[[x_0 v_0],`
\n",
" `[x_1 v_1],`
\n",
" `[x_2 v_2], ...]`
\n",
" To get single list for `x` and single list for `v` we need the transpose of the returned array.
\n",
" (Could also keep return as a single array if that's more useful down the road.)"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [],
"source": [
"t = sp.linspace(0,20,201) # NOTE: The points selected for plotting are \n",
" # not the points used for the numerical \n",
" # evalution.\n",
"x, v = sp.integrate.odeint(damped_osc,u0,t,args=(b,)).T"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [
{
"data": {
"application/javascript": [
"/* Put everything inside the global mpl namespace */\n",
"window.mpl = {};\n",
"\n",
"\n",
"mpl.get_websocket_type = function() {\n",
" if (typeof(WebSocket) !== 'undefined') {\n",
" return WebSocket;\n",
" } else if (typeof(MozWebSocket) !== 'undefined') {\n",
" return MozWebSocket;\n",
" } else {\n",
" alert('Your browser does not have WebSocket support.' +\n",
" 'Please try Chrome, Safari or Firefox ≥ 6. ' +\n",
" 'Firefox 4 and 5 are also supported but you ' +\n",
" 'have to enable WebSockets in about:config.');\n",
" };\n",
"}\n",
"\n",
"mpl.figure = function(figure_id, websocket, ondownload, parent_element) {\n",
" this.id = figure_id;\n",
"\n",
" this.ws = websocket;\n",
"\n",
" this.supports_binary = (this.ws.binaryType != undefined);\n",
"\n",
" if (!this.supports_binary) {\n",
" var warnings = document.getElementById(\"mpl-warnings\");\n",
" if (warnings) {\n",
" warnings.style.display = 'block';\n",
" warnings.textContent = (\n",
" \"This browser does not support binary websocket messages. \" +\n",
" \"Performance may be slow.\");\n",
" }\n",
" }\n",
"\n",
" this.imageObj = new Image();\n",
"\n",
" this.context = undefined;\n",
" this.message = undefined;\n",
" this.canvas = undefined;\n",
" this.rubberband_canvas = undefined;\n",
" this.rubberband_context = undefined;\n",
" this.format_dropdown = undefined;\n",
"\n",
" this.image_mode = 'full';\n",
"\n",
" this.root = $('