{ "metadata": { "name": "" }, "nbformat": 3, "nbformat_minor": 0, "worksheets": [ { "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# PHYS201 PHYSICS IIA - *Python Lab 2*\n", "\n", "*Alexei Gilchrist*, 2014" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# 1 Introduction \n", "\n", "---" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In this lab we will focus on using Python together will some modules to implement matrices and vectors and vartious tasks. We'll start by looking ar the basic object that defines an array in some depth as this is then used throughout the rest of the material." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Following on from the last lab, we will set up the notebook so that we can create plots inline" ] }, { "cell_type": "code", "collapsed": false, "input": [ "%matplotlib inline" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "and import the `matplotlib` routines for plotting" ] }, { "cell_type": "code", "collapsed": false, "input": [ "import matplotlib.pyplot as plt" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Note that here we are using yet another form of the `import` statement. Using the \"`as plt`\" means we can refer to the module by the alias `plt` instead of the lengthy `matplotlib.pyplot`." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# 2 Vectors and Matrices (arrays)\n", "\n", "---" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The nearest structure to a vector that Python has is a list. It really isn't a vector though since vectors are required to have certain properties. To qualify as a vector, or more accurately as an alement of a vector space, a vector has to support two binary operations: addition and multiplication by a scalar. If $\\mathbf{u}$, $\\mathbf{v}$, and $\\mathbf{w}$ are vectors beloging to a vector space and $a$ and $b$ are scalars then the following conditions must be met:\n", "\n", "* Associativity of vector addition: $\\mathbf{u} + (\\mathbf{v}+\\mathbf{w}) = (\\mathbf{u} + \\mathbf{v})+\\mathbf{w}$\n", "* Commutativity: $\\mathbf{u} + \\mathbf{v} = \\mathbf{v}+\\mathbf{u}$\n", "* There exists an identity vector $\\mathbf{0}$ such that $\\mathbf{u} + \\mathbf{0} = \\mathbf{u}$ for all $\\mathbf{u}$\n", "* There exists an inverse $-\\mathbf{u}$ for all $\\mathbf{u}$ such that $\\mathbf{u} + (-\\mathbf{u}) = \\mathbf{0} $\n", "* Associativity of scalar multiplication: $a(b\\mathbf{u}) = (ab)\\mathbf{u}$\n", "* Distributivity of scalar sums: $(a+b)\\mathbf{u} = a\\mathbf{u}+b\\mathbf{u}$\n", "* Distributivity of vector sums: $a(\\mathbf{u}+\\mathbf{v}) = a\\mathbf{u}+a\\mathbf{v}$\n", "* There exists an identity scalar $1$ such that $1\\mathbf{u} = \\mathbf{u}$ for all $\\mathbf{u}$\n", "\n", "The way lists are set up in Python fails these conditions as shown below" ] }, { "cell_type": "code", "collapsed": false, "input": [ "print [1,1,1]+[2,3,4]\n", "print [2,3,4]+[1,1,1]" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "code", "collapsed": false, "input": [ "print 2*([1,2,3]+[4,5,6])\n", "print 2*[1,2,3]+2*[4,5,6]" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "While we *could* create a List object that behaves as required, it's much easier to use someone else's work and simply import some modules! [Numpy](http://numpy.org) contains a multidimensional array object and many routines for doing fast vector and matrix operations. These operations are coded in C and Fotran and are highly optimised making them ideal for numerical computation. Suitably written, 'vectorised' code runs much faster than would be the case in an interpreted language like Python and is comparable to specialised programs like Matlab. \n", "\n", "On top of `numpy` there is another module called [Scipy](http://scipy) which contains additional routines for optimization, special functions, and so on. We'll start by looking at the `array` object from `numpy` then look at other modules later. \n", "\n", "The convention to import `numpy` is to give the module the alias `np`:" ] }, { "cell_type": "code", "collapsed": false, "input": [ "import numpy as np" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 2.1 Vectors" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The `array` object which is like a Python list but contains all the same type of object (e.g. all floating point numbers, or all integers) and can have any number of dimensions. If we create a 1-D array it behaves as a vector. To create it, we pass it a list of values" ] }, { "cell_type": "code", "collapsed": false, "input": [ "u = np.array([1,2,3])" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "code", "collapsed": false, "input": [ "u+u" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "code", "collapsed": false, "input": [ "u-10*u" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can still access and set elements as with standard lists:" ] }, { "cell_type": "code", "collapsed": false, "input": [ "u[2]" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "code", "collapsed": false, "input": [ "u[1]=100\n", "u" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The `array` object will contain the most general data type that you include in the list when creating it. So, for instance, if you pass a list of integers then the array is composed of intergers as in the example above. If you pass a list of floats then you get an array of floats:" ] }, { "cell_type": "code", "collapsed": false, "input": [ "np.array([1.2,3.0,4.0])" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "If there is at least one float and the rest are integers all the integers get converted to floats since in a sense the floats include the integers:" ] }, { "cell_type": "code", "collapsed": false, "input": [ "np.array([1.0,2,3,4,5,6])" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "For the same reason, if you pass a complex number then floats and integers get converted to complex numbers:" ] }, { "cell_type": "code", "collapsed": false, "input": [ "np.array([1j, 2.0, 3, 4,5])" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "You can also pass a second argument to `array` to force the type of object 'l' for long int, 'd' for double float (these are the basic Python number types)" ] }, { "cell_type": "code", "collapsed": false, "input": [ "np.array([2.3,3.,3.4],'l')" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "heading", "level": 3, "metadata": {}, "source": [ " Exercise" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Show each of the coditions for a vector space are met by the `array` object" ] }, { "cell_type": "code", "collapsed": false, "input": [], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 2.2 Matrices" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can use `array` to create matrices too by passing it a list of lists" ] }, { "cell_type": "code", "collapsed": false, "input": [ "M = np.array([[0,1],[1,0]])\n", "u = np.array([3,4])" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "code", "collapsed": false, "input": [ "print M" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Arrays are not set up to do matrix multiplication by default. The method `np.dot` will do matrix-matrix, matrix-vector, and vector-vector multiplication depending on what its arguments are and you will need to use this to do matrix products. There *is* a more specialised array object, `np.matrix`, which will do standard matrix multiplication but it's actually more awkward to use because you have to distinguish column vectors from row vectors, and many functions return `array`s by default. `array` is more general too as it allows any number of dimensions so we'll stick to just using `array` objects." ] }, { "cell_type": "code", "collapsed": false, "input": [ "np.dot(M,u)" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "heading", "level": 3, "metadata": {}, "source": [ " Exercise" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Create a 3-D array" ] }, { "cell_type": "code", "collapsed": false, "input": [], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 2.3 Array methods" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "An `array` object has many useful methods" ] }, { "cell_type": "code", "collapsed": false, "input": [ "M = np.array([[1,2,3],[4,5,6],[7,8,9]]) # set up a 2-D array for demonstration\n", "print(M)" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "####Transpose" ] }, { "cell_type": "code", "collapsed": false, "input": [ "print M.transpose()\n", "print\n", "print M.T # this is a shortcut for .transpose()" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "####Trace" ] }, { "cell_type": "code", "collapsed": false, "input": [ "print M.trace() # the sum of the diagonal" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Diagonal of matrix" ] }, { "cell_type": "code", "collapsed": false, "input": [ "print M.diagonal()" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### size and shape" ] }, { "cell_type": "code", "collapsed": false, "input": [ "print M.size\n", "print M.shape" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### min, max, mean, sum, prod (product), var (variance), std (standard deviation)" ] }, { "cell_type": "code", "collapsed": false, "input": [ "print \"1. min: \", M.min()\n", "print \"2. max: \", M.max()\n", "print \"3. mean:\", M.mean()\n", "print \"4. sum: \", M.sum()\n", "print \"5. prod:\", M.prod()\n", "print \"6. std: \", M.std()\n", "print \"7. var: \", M.var()" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "####Conjugate, real and imaginary parts" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "First create a complex array from the one we have" ] }, { "cell_type": "code", "collapsed": false, "input": [ "M2 = M*1j\n", "print M2" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "code", "collapsed": false, "input": [ "print M2.conj()\n", "print\n", "print M2.real\n", "print\n", "print M2.imag" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "There are also other methods that you can look up" ] }, { "cell_type": "heading", "level": 3, "metadata": {}, "source": [ " Exercise" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Show how to perform the conjugate transpose operation by chaining some of the operations above" ] }, { "cell_type": "code", "collapsed": false, "input": [], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 2.4 Creating arrays" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The following two methods are pretty self evident, the argument is the shape or size of the array to create" ] }, { "cell_type": "code", "collapsed": false, "input": [ "np.ones([2,3]) # all 1's" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "code", "collapsed": false, "input": [ "np.zeros([3,2]) # all 0's" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "code", "collapsed": false, "input": [ "np.eye(4) # the identity" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "To create a sequence of numbers use `arange` (i.e. array range). This works the same as the function we created in the first lab, only it's much more efficient and returns an array instead of a Python list." ] }, { "cell_type": "code", "collapsed": false, "input": [ "np.arange(4,20,2) # (start, end, step-size) end point is not inclusive" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Alternatively you can use `linspace` to specify a start, and end and the number of points you want. " ] }, { "cell_type": "code", "collapsed": false, "input": [ "np.linspace(0,10,21) # (start, end, number of points) end point is inclusive" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "(think about why it's 21 points to get the spacing of 0.5 for a moment)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can also create an array with a particular sequence on the diagonal with `diag`:" ] }, { "cell_type": "code", "collapsed": false, "input": [ "np.diag(np.arange(4))" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "heading", "level": 3, "metadata": {}, "source": [ " Exercise" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Look up the documentation for np.tile (using either `help` or putting a \"?\" after the command) and use it to create an array like the one below:\n", "\n", " [[1, 2, 1, 2],\n", " [3, 4, 3, 4],\n", " [1, 2, 1, 2],\n", " [3, 4, 3, 4]]" ] }, { "cell_type": "code", "collapsed": false, "input": [], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 2.5 Functions on arrays" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "`numpy` also defines many common mathematical functions which will efficiently evaluate on the entire array. For example, we can generate a range in a variable $x$ and then plot $\\sin(x)$." ] }, { "cell_type": "code", "collapsed": false, "input": [ "x = np.linspace(0,2*np.pi,50)\n", "plt.plot(x,np.sin(x))" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "`numpy` has automatically looped the function `sin( )` over the array `x` and returned a new array of the results." ] }, { "cell_type": "heading", "level": 3, "metadata": {}, "source": [ " Exercise" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Construct an array of the cubes of the first 100 integers" ] }, { "cell_type": "code", "collapsed": false, "input": [], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 2.6 Array conditions and tests" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Comparisons with arrays and a value are vecorised across the array and return an array of booleans" ] }, { "cell_type": "code", "collapsed": false, "input": [ "u = np.arange(10) # set up an array for examples\n", "print u" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "code", "collapsed": false, "input": [ "u == 6" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "code", "collapsed": false, "input": [ "u > 7" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "code", "collapsed": false, "input": [ "u**2 > 20" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 2.7 Indexing" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The usual Python list indexing notation are available for arrays (the notation is `start:end:step`)" ] }, { "cell_type": "code", "collapsed": false, "input": [ "u = np.arange(10,20) # set up an array for examples\n", "print u" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "code", "collapsed": false, "input": [ "print u[0] # first element (remember, indexes start with 0)\n", "print u[1:3] \n", "print u[-1] # last element\n", "print u[::2] # every 2nd element" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "These conventions are extended to multi-dimentional arrays by having multiple arguments to the index. Here are a bunch of examples, make sure you understand how they work" ] }, { "cell_type": "code", "collapsed": false, "input": [ "M = np.array([[1,2,3,4],[5,6,7,8],[9,10,11,12],[13,14,15,16]]) # set up array for examples\n", "print M" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "code", "collapsed": false, "input": [ "M[1,1] # a single element" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "code", "collapsed": false, "input": [ "M[:,2] # a column" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "code", "collapsed": false, "input": [ "M[0:2,:] # some rows" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "code", "collapsed": false, "input": [ "M[1:,1:] # a sub-array" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Something to be aware of is that in `numpy`, indexing creates a *view* into the array not a copy. So if you modify the view, you are modifying the original at the same time. To get a copy of the array use the `copy()` method. The example below demonstrates " ] }, { "cell_type": "code", "collapsed": false, "input": [ "M2 = M.copy()\n", "N = M2[:,2]\n", "print \"N=\", N\n", "N[0] = 0\n", "print \"N=\", N\n", "print \"M2=\", M2\n", "print \"M=\", M" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "A fancier form of indexing is to use an array of booleans as a mask, the result will be an array (a copy not a view) with the items at the locations where the indexing array is `True`. This is especially powerful as we can generate an array of booleans with a test." ] }, { "cell_type": "code", "collapsed": false, "input": [ "# set up an arrays for examples\n", "u = 2*np.arange(10)\n", "div3 = u%3==0 # this will be True when item in u is divisible by 3\n", "print u\n", "print div3" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "code", "collapsed": false, "input": [ "u[div3] # this will pick out the items divisible by 3" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "A compact way of doing the above is to generate the array of booleans from the test on the fly:" ] }, { "cell_type": "code", "collapsed": false, "input": [ "u[u%3==0]" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Finally, you can also pass an array a list of indexes to pick out. These can be repeated and in any order:" ] }, { "cell_type": "code", "collapsed": false, "input": [ "u[[1,1,9,1]]" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "As you can see, `numpy` is very flexible with how to index an array. This provides a surprising amount of power as often quite complicated tasks can be solved by a clever use of indexing. " ] }, { "cell_type": "heading", "level": 3, "metadata": {}, "source": [ " Exercise" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Create an array of the first 100 numbers. By using indexing, set to zero all those numbers not divisible by 3." ] }, { "cell_type": "code", "collapsed": false, "input": [], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# 3 Linear Algebra" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "`numpy` has a linear algebra sub module that has methods for computing determinants, eigenvalues etc. For convenience we'll import it and give it the alias `la`:" ] }, { "cell_type": "code", "collapsed": false, "input": [ "import numpy.linalg as la" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Check the documentation to see what is defined with this handy trick (you can close the frame with the X on the right):" ] }, { "cell_type": "code", "collapsed": false, "input": [ "la?" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In this lab we will only make use of the eigenvalue methods. The method `la.eig(M)` will compute the eigenvalues and eigenvectors of a matrix `M` (a 2-D array). The eigenvalues are returned in an array followed by the eigenvectors as the columns of a 2-D array." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In the solved problems for the unit, there is a coupled oscillator problem which has two masses suspended in a chain by springs. The top mass is $m$ and the bottom mass is $2m$. Both springs have spring constant $k$. In the solution the system can be reduced to the matrix equation\n", "\n", "\\begin{equation}\n", "\\frac{k}{m}\n", "\\begin{pmatrix}\n", "2 & -1 \\\\\n", "-1/2 & 1/2\n", "\\end{pmatrix} \n", "\\begin{pmatrix}\n", "A \\\\ B\n", "\\end{pmatrix}= \\omega^2\n", "\\begin{pmatrix}\n", "A \\\\ B\n", "\\end{pmatrix}\n", "\\end{equation}\n", "\n", "and the eigenvalues are calculated to be $(5k\\pm\\sqrt{17}k)/(4m)$ and the eigenvectors are $\\begin{pmatrix}\n", "(-3\\pm\\sqrt{17})/2 \\\\ 1\n", "\\end{pmatrix}$." ] }, { "cell_type": "heading", "level": 3, "metadata": {}, "source": [ " Exercise" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Verify this solution using `la.eig`. N.B. use floating point numbers when you write divisors (e.g. `1/2.0`) or Python's integer division will bite you! Normalising the vectors given in the example above will help you to verify the solution." ] }, { "cell_type": "code", "collapsed": false, "input": [], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# 4 Dynamical system maps" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "we can reduce an $n^\\mathrm{th}$ order differential equation into $n$ first order differential equations:\n", "\n", "\\begin{align}\n", " \\dot{x}_1 &= f_1(x_1, x_2, \\ldots) \\\\\n", " \\dot{x}_2 &= f_2(x_1, x_2, \\ldots) \\\\\n", " & \\vdots\n", "\\end{align}\n", "\n", "These equations can be viewed as giving the components of an $n$-dimensional arrow pointing in the direction of change. Where the arrows are length zero defines *fixed points*. Finding these is a matter of solving the set of equations obtained by setting $\\dot{x}_1 = \\dot{x}_2 = \\ldots = 0$. But what about plotting the direction of the arrows? This is tedious by hand, but fortunately easy to do for two dimensions with `plplot.quiver`." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "A rigid-rod pendulum has the equation of motion\n", "\n", "\\begin{equation}\n", " \\ddot{\\theta} + \\omega_0^2 \\sin(\\theta) = 0\n", "\\end{equation}\n", "\n", "This can be trivially converted to two first order differential equations:\n", "\n", "\\begin{align}\n", " \\dot{\\theta} &= v \\\\\n", " \\dot{v} &= -\\omega_0^2 \\sin(\\theta)\n", "\\end{align}\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "To plot a grid of arrows $(\\dot{\\theta},\\dot{v})$ we create arrays for $\\theta$ and $v$ that span the range we are interested in, convert them to the matrices that are needed for the plot routine. We can then efficiently compute $\\dot{\\theta}$ and $\\dot{v}$ using `numpy`'s element wise operations. See the example below." ] }, { "cell_type": "code", "collapsed": false, "input": [ "w0=1.0\n", "\n", "# create regularly spaced arrays in x and v\n", "x = np.linspace(-2*np.pi,2*np.pi,15) \n", "v = np.linspace(-2,2,15)\n", "\n", "# convert them to matrices with repeated rows for X and repeated columns for V\n", "# (this format is needed by the plt.quiver plot)\n", "X, V = np.meshgrid(x, v)\n", "\n", "# calculate the derivatives for x and v - these are the components of the arrows\n", "# Note that numpy will do the calculation for each element of the arrays\n", "DX = V\n", "DV = -w0**2*np.sin(X)" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now we can plot the arrows with `plplot.quiver`" ] }, { "cell_type": "code", "collapsed": false, "input": [ "plt.figure(figsize=(8, 4)) # create a figure \n", "plt.quiver(X,V,DX,DV,alpha=.5)" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "heading", "level": 3, "metadata": {}, "source": [ " Exercise" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Add damping of $\\gamma=1.5$ to the above model and plot the dynamical map." ] }, { "cell_type": "code", "collapsed": false, "input": [], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "heading", "level": 3, "metadata": {}, "source": [ " Exercise" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Plot the dynamical map for the nonlinear oscillator in the range $-2\\le x \\le 2$ and $-2\\le v \\le 2$.\n", "\n", "\\begin{equation}\n", " \\ddot{x} + b x^3 - c x = 0\n", "\\end{equation}\n", "\n", "for the parameter $b=1$ and $c=2$." ] }, { "cell_type": "code", "collapsed": false, "input": [], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# 5 Numerical Integration" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "the module `scipy.integrate` provides a bunch of functions for doing numerical integration. We'll use the function `odeint` from that module." ] }, { "cell_type": "code", "collapsed": false, "input": [ "from scipy.integrate import odeint" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Again, we can reduce an $n^\\mathrm{th}$ order differential equation into $n$ first order differential equations which we can write compactly as $\\dot{\\mathbf{x}} = f(\\mathbf{x},t)$, where $\\mathbf{x}$ is a vector of variables and $f$ is some function. We can numerically solve such functions by using `odeint`:\n", "\n", " xt = odeint(f, x0, t, args=(...))\n", " \n", "Here, `f` is a python function that depends on an array of values and time: `f(x,t, ...)` and it returns an array of the time derivatives, `x0` is an array of initial conditions, `t` is an array of the times at which you want the solution, and `xt` is a 2-D array of the results of `x` at those times. Note that the objects in the keyword argument `args=` in `odeint` are passed to the method `f` after the `x` and `t` as indicated. This gives us the ability to pass parameters into the function, it can be omitted if the function doesn't take any extra parameters." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "For example, the damped harmonic oscillator has the equation of motion\n", "\n", "\\begin{equation}\n", "\\ddot{x}+\\gamma \\dot{x}+\\omega_0^2 x = 0\n", "\\end{equation}\n", "\n", "Which can we written as two first order equations:\n", "\n", "\\begin{align}\n", " \\dot{x} &= v \\\\\n", " \\dot{v} &= -\\gamma v-\\omega_0^2 x\n", "\\end{align}\n", "\n", "Or in vector notation as $\\dot{\\mathbf{x}} = f(\\mathbf{x},t)$ where $\\mathbf{x}=\\begin{pmatrix}x\\\\ v\\end{pmatrix}$, and $f(\\mathbf{x},t)\\equiv \\begin{pmatrix}\\dot{x}\\\\ \\dot{v}\\end{pmatrix} = \\begin{pmatrix}v \\\\ -\\gamma v+\\omega_0^2 x\\end{pmatrix}$ ." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "First define the function $f$:" ] }, { "cell_type": "code", "collapsed": false, "input": [ "def fdampedHO(x, t, gamma, omega0):\n", " derivatives = np.array([ x[1], -gamma*x[1]-(omega0**2)*x[0] ])\n", " return derivatives" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Then the initial conditions:" ] }, { "cell_type": "code", "collapsed": false, "input": [ "x0 = np.array([1,0])\n", "t = np.linspace(0,30,200)" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now solve the system for various parameter regimes:" ] }, { "cell_type": "code", "collapsed": false, "input": [ "xt1 = odeint(fdampedHO, x0, t, args=(0,1)) # undamped\n", "xt2 = odeint(fdampedHO, x0, t, args=(0.5,1)) # under-damped\n", "xt3 = odeint(fdampedHO, x0, t, args=(2,1)) # critically damped\n", "xt4 = odeint(fdampedHO, x0, t, args=(4,1)) # over-damped" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Plot everything. Note that we are interested in only the position component of `xt[:,0]`, (the other component `xt[:,1]` is the velocity)." ] }, { "cell_type": "code", "collapsed": false, "input": [ "fig, ax = plt.subplots()\n", "ax.plot(t, xt1[:,0], 'k', label=\"undamped\", linewidth=0.25)\n", "ax.plot(t, xt2[:,0], 'r', label=\"under-damped\")\n", "ax.plot(t, xt3[:,0], 'b', label=\"critically damped\")\n", "ax.plot(t, xt4[:,0], 'g', label=\"over-damped\")\n", "ax.legend();" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We could plot the position against the velocity too" ] }, { "cell_type": "code", "collapsed": false, "input": [ "fig, ax = plt.subplots()\n", "ax.plot(xt1[:,0], xt1[:,1], 'k', label=\"undamped\", linewidth=0.25)\n", "ax.plot(xt2[:,0], xt2[:,1], 'r', label=\"under-damped\")\n", "ax.plot(xt3[:,0], xt3[:,1], 'b', label=\"critically damped\")\n", "ax.plot(xt4[:,0], xt4[:,1], 'g', label=\"over-damped\")\n", "ax.legend();" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "heading", "level": 3, "metadata": {}, "source": [ " Exercise" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Modify the above code to calculate and redo the plots for the *driven* and damped oscillator:\n", "\n", "\\begin{align}\n", " \\dot{x} &= v \\\\\n", " \\dot{v} &= -\\gamma v-\\omega_0^2 x + f_0 \\cos(\\omega_f t)\n", "\\end{align}\n", "\n", "for the parameters $\\gamma = 0.5$, $\\omega_0=2.0$, $\\omega_f=2.3$, and $f_0=0.1$." ] }, { "cell_type": "code", "collapsed": false, "input": [], "language": "python", "metadata": {}, "outputs": [] } ], "metadata": {} } ] }