{ "metadata": { "name": "" }, "nbformat": 3, "nbformat_minor": 0, "worksheets": [ { "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# PHYS201 PHYSICS IIA - *Python Lab 3*\n", "\n", "*Alexei Gilchrist*, 2014\n", "\n", "---" ] }, { "cell_type": "heading", "level": 1, "metadata": {}, "source": [ "1 Introduction " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In this lab we are going to use the Python module sympy to do *symbolic* computation. This is quite different from the last lab were we where principally dealing with arrays of numbers and manipulating them. In this lab we will keep things as symbols and manipulate the symbols to perform algebra and calculus. This is very similar to commercial programs such as Mathematica.\n", "\n", "The sympy module is a pure-python module and is not really optimised for big calculations. It is sufficiently powerful to handle most of the problems encountered in phys201 though, and it's an alternative to using Mathematica or [Wolfram Alpha](http://www.wolframalpha.com). We won't touch some of the more specialised sub modules, just sticking to the basics. Think of it as a super-charged calculator." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Some resources for sympy:\n", "\n", "* [Main documentation](http://docs.sympy.org/latest/index.html) - official documentation\n", "* [Live Sympy shell](http://live.sympy.org) - run code here if you need to do a quick calculation\n", "* [SymPy Tutorial](http://docs.sympy.org/latest/tutorial/index.html) - this tutorial has live code-blocks you can modify.\n", "* [Sympy - Symbolic algebra in Python](http://nbviewer.ipython.org/github/jrjohansson/scientific-python-lectures/blob/master/Lecture-5-Sympy.ipynb) - tutorial by J.R. Johansson\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In this lab we'll also take a different tack from the previous labs in that we will import everything in the module to the local namespace with from sympy import *. This will provide a more convenient calculation environment but keep in mind that many of the methods come from sympy not from native Python." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "*Remember to execute the cells as you go!*" ] }, { "cell_type": "code", "collapsed": false, "input": [ "from sympy import *" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "One of the great things about sympy is that the output can be formatted with the typesetting language $\\LaTeX$. To have this on by default run the following:" ] }, { "cell_type": "code", "collapsed": false, "input": [ "init_printing(use_latex=True)" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now, for example, the following displays nicely" ] }, { "cell_type": "code", "collapsed": false, "input": [ "(1+pi)**2" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "sympy will keep things as symbols by default, so that while the following should not be surprising:" ] }, { "cell_type": "code", "collapsed": false, "input": [ "sqrt(9)" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "but this is pretty neat:" ] }, { "cell_type": "code", "collapsed": false, "input": [ "sqrt(8)" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "heading", "level": 1, "metadata": {}, "source": [ "2 Basics" ] }, { "cell_type": "heading", "level": 2, "metadata": {}, "source": [ "2.1 Symbols" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "To create objects that behave as algebraic symbols use the Symbol class. " ] }, { "cell_type": "code", "collapsed": false, "input": [ "x = Symbol('x')" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "or the shortcut method symbols which can return several symbol objects in one go:" ] }, { "cell_type": "code", "collapsed": false, "input": [ "x, y, z, t = symbols('x y z t')" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "These objects now behave as you'd expect from other symbolic computation programs:" ] }, { "cell_type": "code", "collapsed": false, "input": [ "x+3*y-(2+x)**2" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Generally it is really a good idea to keep the variable name the same as the symbol itself like was done above. But there are situations when you might not want to do this, for intance we can give the symbol a nice representation in $\\LaTeX$:" ] }, { "cell_type": "code", "collapsed": false, "input": [ "omega0 = Symbol('\\omega_0')" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "code", "collapsed": false, "input": [ "omega0**2" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "As another shortcut, we can import pre-made symbols from sympy.abc which defines the following" ] }, { "cell_type": "code", "collapsed": false, "input": [ "import sympy.abc\n", "dir(sympy.abc)" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "code", "collapsed": false, "input": [ "from sympy.abc import alpha, w" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "code", "collapsed": false, "input": [ "w+alpha" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Keep in mind that these sympols are special objects. So if you run something like x=4, x will no longer by a symbolic variable it will just be the object '4'. Also python will not prevent you from overriding default sympy names or functions (which we imported with from sympy import *), so be careful. As a rule it is recommended that you not use I, E, S, N, C, O, or Q for variable or symbol names as those are used by sympy." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can however create a variable which holds a symbolic expression using the normal Python syntax" ] }, { "cell_type": "code", "collapsed": false, "input": [ "z2 = omega0+3" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "code", "collapsed": false, "input": [ "z2**2" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Symbols can also have certain assumptions attached to them like that they are real, or positive, or integers etc" ] }, { "cell_type": "code", "collapsed": false, "input": [ "n = Symbol('n', positive=True)" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "code", "collapsed": false, "input": [ "n > 0" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "heading", "level": 3, "metadata": {}, "source": [ " Exercise" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Just to be perverse, define a symbol a that displays as b, and a symbol b that displays a and do a computation with them." ] }, { "cell_type": "code", "collapsed": false, "input": [], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "heading", "level": 2, "metadata": {}, "source": [ "2.2 Numbers" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "sympy has a different set of objects to represent numbers than Python does. It uses an artitrary precision library for computation and defines mathematical constants such as pi, e, and oo (infinity). It uses I$=\\sqrt{-1}$, Real, Rational, and Integer to represent numbers, of which the two to be aware of is I and Rational as the other two stay mostly behind the scenes." ] }, { "cell_type": "code", "collapsed": false, "input": [ "1/oo" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "heading", "level": 3, "metadata": {}, "source": [ "2.2.1 Complex numbers" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "sympy uses I for a complex number" ] }, { "cell_type": "code", "collapsed": false, "input": [ "3*(2+I)" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "heading", "level": 3, "metadata": {}, "source": [ "2.2.2 Rational numbers" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "By defining rational numbers as special objects, we can manipulate them algebraically instead of reducing them to floating point numbers." ] }, { "cell_type": "code", "collapsed": false, "input": [ "r1 = Rational(4,5)\n", "r2 = Rational(5,4)\n", "print r1, r2" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "code", "collapsed": false, "input": [ "r1/r2" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "code", "collapsed": false, "input": [ "r1+r2" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Compare the following results" ] }, { "cell_type": "code", "collapsed": false, "input": [ "int(2)/int(3)" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "code", "collapsed": false, "input": [ "Integer(2)/Integer(3)" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "heading", "level": 3, "metadata": {}, "source": [ " Exercise" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Calculate $i^{1/3}$ symbolically " ] }, { "cell_type": "code", "collapsed": false, "input": [], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "heading", "level": 3, "metadata": {}, "source": [ "2.2.3 Numerical evaluation" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "evalf or N can be used to evaluate an expression or constant explicitly. the method can take an optional argument for the number of significant digits" ] }, { "cell_type": "code", "collapsed": false, "input": [ "print sqrt(8)\n", "print sqrt(8).evalf()\n", "print N(sqrt(8))" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "heading", "level": 3, "metadata": {}, "source": [ " Exercise" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Look up the help for N and print off the first thousand digits of $\\pi$" ] }, { "cell_type": "code", "collapsed": false, "input": [], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "heading", "level": 2, "metadata": {}, "source": [ "2.3 Equals" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Sometimes, seemingly simple things can turn around and bite you. What do we mean when we write an expression like\n", "$(x-1)^2 = x^2 - 2x+1$? We've already come across two uses of the symbol '=' in Python. A single equals (=) assigns some object to a variable, and a double equals (==) compares two objects following some fairly complicated [set of rules](https://docs.python.org/2/reference/expressions.html#not-in). Neither of these uses are what we want in algebra or calculus, what we want is to say that the expression on the left is the same as the expression on the right subject to some manipulation. " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now the = symbol is hardwired into python so we can't do much with that, and == has a specific meaning for sympy expressions - it means that the expressions and the left and right are *structurally* the same:" ] }, { "cell_type": "code", "collapsed": false, "input": [ "print (x-1)**2 == x**2-2*x+1 # structures differ\n", "print x**2-2*x+1 == x**2-2*x+1 # structures are the same" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "So where does this leave us? Well, there are two functions we might want from an = in a computer algebra system. One is to make a statement that two expressions should be taken to be equal like for an equation. For this we can use Eq (which is a synonym for Equality):" ] }, { "cell_type": "code", "collapsed": false, "input": [ "Eq(x**2-2*x+1,0)" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The other sense is that we want to *test* if two expressions are the same. Surprisingly, it turns out this is not possible in general! (see [Richardson's theorem](http://en.wikipedia.org/wiki/Richardson%27s_theorem)). What we can do is try simplify both sides into a particular form, or turn $A=B$ into $A-B=0$ and try and reduce $A-B$ to 0. For this task sympy has some methods like simplify and expand that we can use." ] }, { "cell_type": "code", "collapsed": false, "input": [ "expand((x-1)**2) == x**2-2*x+1" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "code", "collapsed": false, "input": [ "simplify( (x**2-2*x+1) - (x-1)**2 )" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "There is also equals method which tests the expressions at random points" ] }, { "cell_type": "code", "collapsed": false, "input": [ "( (x-1)**2 ).equals( x**2-2*x+1 )" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "heading", "level": 2, "metadata": {}, "source": [ "2.4 Algebra" ] }, { "cell_type": "heading", "level": 3, "metadata": {}, "source": [ "2.4.1 Substitution" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "If we have some expression" ] }, { "cell_type": "code", "collapsed": false, "input": [ "expr = x**2-2*x+1\n", "expr" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "values can be substituted into it using subs" ] }, { "cell_type": "code", "collapsed": false, "input": [ "expr.subs(x,3)" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Pass a list of tuples to it to expand multiple symbols" ] }, { "cell_type": "code", "collapsed": false, "input": [ "(x**2+y**3).subs([(x,1),(y,3)])" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "or a dict works just as well" ] }, { "cell_type": "code", "collapsed": false, "input": [ "(x**2+y**3).subs({x:1, y:3})" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Substitution is not limited to numerical values! By substituting in varios expressions you can perform algebraic manipulations" ] }, { "cell_type": "code", "collapsed": false, "input": [ "expr.subs(x, x+1)" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "heading", "level": 3, "metadata": {}, "source": [ "2.4.2 Expand and factor" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "expand multiplies out factors and factor does factoring:" ] }, { "cell_type": "code", "collapsed": false, "input": [ "expr = (x+1)*(x+3)*(x-4)\n", "expr" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "code", "collapsed": false, "input": [ "expr2=expand(expr)\n", "expr2" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "code", "collapsed": false, "input": [ "factor(expr2)" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "It's often possible to pass hints to functions to tell them what it is you want to do, for instance trig=True will tell expand to try and expand trigonometric functions. See expand?, factor? to get more information." ] }, { "cell_type": "code", "collapsed": false, "input": [ "print sin(x+y), \"=\", expand( sin(x+y) , trig=True)\n", "print sin(x-y), \"=\", expand( sin(x+y) , trig=True)" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "heading", "level": 3, "metadata": {}, "source": [ "2.4.3 Simplify" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "One of the useful tasks of a computer algebra system is to simplify expressions. There are many of functions to perform this task in sympy, and one general function called simplify that attempts to cleverly apply what is necessary." ] }, { "cell_type": "code", "collapsed": false, "input": [ "simplify(sin(x)**2 + cos(x)**2)" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "code", "collapsed": false, "input": [ "simplify(cos(x)/sin(x))" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "There are also more specialised finctions you can use: trigsimp, powsimp, logcombine, etc. See [Simplification](http://docs.sympy.org/latest/tutorial/simplification.html) for more examples and techniques." ] }, { "cell_type": "heading", "level": 3, "metadata": {}, "source": [ "2.4.4 Apart and together" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Expressions with fractions can be manipulated with apart and together:" ] }, { "cell_type": "code", "collapsed": false, "input": [ "expr = 1/(x+2)+1/(x+3)" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "code", "collapsed": false, "input": [ "expr2=together(expr)\n", "expr2" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "code", "collapsed": false, "input": [ "apart(expr2)" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "heading", "level": 2, "metadata": {}, "source": [ "2.5 Calculus" ] }, { "cell_type": "heading", "level": 3, "metadata": {}, "source": [ "2.5.1 Differentiation" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Differentiation can be performed by diff. For example, if we set $y$ to the solution for a critically damped harmonic oscillator," ] }, { "cell_type": "code", "collapsed": false, "input": [ "from sympy.abc import gamma, A, B # set some symbols we are going to use\n", "y = (A+B*t)*exp(-gamma/2*t)\n", "y" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "It's easy to find the velocity by differentiating (and simplifying)," ] }, { "cell_type": "code", "collapsed": false, "input": [ "simplify(diff(y,t))" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "and the acceleration by differentiating twice," ] }, { "cell_type": "code", "collapsed": false, "input": [ "simplify(diff(y,t,2))" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "heading", "level": 3, "metadata": {}, "source": [ " Exercise" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "A under-damped harmonic oscillator has the solution\n", "\n", "\$$\n", " x(t) = A e^{-\\gamma/2 t}\\cos(\\omega_d t)\n", "\$$\n", "\n", "Differentiate to calculate the velocity and the acceleration. (don't forget to define suitable symbols)" ] }, { "cell_type": "code", "collapsed": false, "input": [], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "heading", "level": 3, "metadata": {}, "source": [ "2.5.2 Integration" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "And of course, there is also integration. You can do indefinite:" ] }, { "cell_type": "code", "collapsed": false, "input": [ "integrate(exp(-x**2),x)" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "definite:" ] }, { "cell_type": "code", "collapsed": false, "input": [ "integrate(exp(-x**2),(x,-oo,oo))" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "and multiple integrals:" ] }, { "cell_type": "code", "collapsed": false, "input": [ "x,y = symbols('x y')\n", "integrate(exp(-x**2-2*y**2),(x,-oo,oo),(y,-oo,oo))" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "heading", "level": 3, "metadata": {}, "source": [ " Exercise" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Imagine rotating a square slab of material about an axis perpendicular to the slab. The slab has mass $M$ side length $L$. The moment of inertial for this system is \n", "\$$\n", "I = \\int r^2 dm = \\frac{M}{L^2} \\iint (x^2+y^2) dx dy\n", "\$$\n", "since for this shape $r=\\sqrt(x^2+y^2)$ and a small element of mass is $dm = \\frac{M}{L^2}dxdy$\n", "\n", "By setting different limits on the integrals, calculate the moment of inertia for an axis that passes through\n", "\n", "1. the middle of the square slab\n", "1. a corner\n", "1. the middle of a side" ] }, { "cell_type": "code", "collapsed": false, "input": [], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "heading", "level": 3, "metadata": {}, "source": [ "2.5.3 Series" ] }, { "cell_type": "code", "collapsed": false, "input": [ "series(sin(x), x)" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "This is equivalent" ] }, { "cell_type": "code", "collapsed": false, "input": [ "sin(x).series(x)" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Or to expand it about $x_0$ and include to 8th order" ] }, { "cell_type": "code", "collapsed": false, "input": [ "x0 = Symbol(\"x_0\")\n", "sin(x).series(x, x0, 8)" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "heading", "level": 3, "metadata": {}, "source": [ " Exercise" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "It's common to make the aproximation $(1+x)^{1/2} \\approx 1+\\frac{1}{2}x$. \n", "\n", "* Use series to show why this approximation works. \n", "* What is $\\frac{1}{\\sqrt{1-v^2}}$ to 4th order?" ] }, { "cell_type": "code", "collapsed": false, "input": [], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "heading", "level": 2, "metadata": {}, "source": [ "2.6 Linear Algebra" ] }, { "cell_type": "heading", "level": 3, "metadata": {}, "source": [ "2.6.1 Matrices" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Matrices can also be constructed in sympy in a similar way to numpy except that they can now have symbols as entries. The basic object is the Matrix which takes a list of lists. There are differences with numpy though which can be annoying." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "For example, here is a matrix that turns up in a particular coupled oscillator problem: " ] }, { "cell_type": "code", "collapsed": false, "input": [ "k,m = symbols('k m')\n", "M = k/m*Matrix([[1, -Rational(1,2)],[-1, 2]])\n", "M" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "There are various ways of creating matrices just as there are in numpy:" ] }, { "cell_type": "code", "collapsed": false, "input": [ "print eye(3) # identity matrix 3x3\n", "print\n", "print zeros(2,3) # full of zeros 2x3\n", "print \n", "print ones(3,1) # column vector of ones\n", "print \n", "print diag(1,2,3) # 1,2,3 on diag (NB this is different from numpy which takes a list!)\n", "print\n", "print diag(ones(2),k*ones(3),2) # sympy's diag allows block diagonal construction" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "And you can take the transpose or get certain columns or vectors with standard notation:" ] }, { "cell_type": "code", "collapsed": false, "input": [ "M2 = Matrix([[1, 2, 3], [4, 5, 6]])\n", "print M2\n", "print\n", "print M2.T # transpose\n", "print\n", "print M2[1,:] # 2nd row\n", "print\n", "print M2[:,0] # first column" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "But of course functions are calculated symbolically which can make sympy much more useful. The following calculate the determinant and the eigenvalues and eigenvectors:" ] }, { "cell_type": "code", "collapsed": false, "input": [ "M.det()" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "code", "collapsed": false, "input": [ "print M.eigenvals()" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "code", "collapsed": false, "input": [ "M.eigenvects()" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Note that the eigenvalues and eigenvectors are a bit buried in the structure that is returned and you have to dig them out to use them further. eigenvects returns a list, each element is also a list made up of (*eigenvalue*, *duplicity*, *eigenvector*). The eigenvector in it's own list too!" ] }, { "cell_type": "heading", "level": 3, "metadata": {}, "source": [ " Exercise" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "If we couple two pendulums by a spring with one of the pendulums three times the mass of the other, the resulting equations of motion of the form\n", "\\begin{align}\n", "\\ddot{x}&+bx +a(x-y) =0 \\\\\n", "\\ddot{y}&+by +3a(y-x) =0.\n", "\\end{align}\n", "\n", "The trial solutions $x=Ae^{i\\omega t}$ and $y=Be^{i\\omega t}$ lead to the following equations\n", "\\begin{align}\n", "bA &+a(A-B) =\\omega^2 A\\\\\n", "bB &+3a(B-A) =\\omega^2 B.\n", "\\end{align}\n", "\n", "* Define $a$ and $b$ as symbols and write the resulting matrix $M$ so that the system can be written as the eigenvalue equation $Mv=\\omega^2v$ with $v=[A,B]^T$. \n", "* Find the eigenvalues and eigenvectors of the matrix.\n", "* set e1, and e2 as the eigenvalues, v1, and v2 as the eigenvectors. (You will need to index into the returned results)" ] }, { "cell_type": "code", "collapsed": false, "input": [], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "heading", "level": 3, "metadata": {}, "source": [ " Exercise" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In the above system of equations the following where defined: $a=k/(3m)$ and $b=g/l$\n", "\n", "* create the symbols $g$, $l$, $k$, and $m$ defining them to be positive (this will help with simplification)\n", "* substitute in the expressions for $a$ and $b$ using subs\n", "* simplify the result if necessary" ] }, { "cell_type": "code", "collapsed": false, "input": [], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "heading", "level": 3, "metadata": {}, "source": [ "2.6.2 Solving Equations" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "One of the most common things to do is solve equations. For example, solving for $x$ in $x^2-1=0$. Although this equation is trivial we'll use sympy to solve it to see how it works. First we'll define some symbols then define an equation" ] }, { "cell_type": "code", "collapsed": false, "input": [ "x,y = symbols('x y') \n", "eq1 = Eq( x**2-1,0)\n", "eq1" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now tell sympy to solve the equation for $x$:" ] }, { "cell_type": "code", "collapsed": false, "input": [ "solve(eq1,x)" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In this intance we can simplify the process and give only the LHS and sympy will assume the RHS is zero:" ] }, { "cell_type": "code", "collapsed": false, "input": [ "solve(x**2-1 ,x)" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "solve can also takle systems of equations such as\n", "\\begin{align}\n", "x+y-1 & = 0 \\\\\n", "x-y-1 & = 0\n", "\\end{align}" ] }, { "cell_type": "code", "collapsed": false, "input": [ "x,y = symbols('x y')\n", "solve([x + y - 1, x - y - 1], [x,y])" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "heading", "level": 3, "metadata": {}, "source": [ " Exercise" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The formula\n", "\$$\n", "\\begin{pmatrix} x \\\\ y \\end{pmatrix} = X v_1 + Y v_2\n", "\$$\n", "defines the physical variables $x$ and $y$ in terms of the normal mode variables $X$ and $Y$ using the eigenvectors.\n", "\n", "* Define symbols X and Y\n", "* Create a pair of equations using Eq for the eigenvectors found in the previous exercise\n", "* solve the equations for X and Y" ] }, { "cell_type": "code", "collapsed": false, "input": [], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "heading", "level": 3, "metadata": {}, "source": [ " Exercise" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The Lotka-Voltera equations are \n", "\\begin{align}\n", "\\dot{x}&= Ax-Bxy \\\\\n", "\\dot{y}&= Dxy-Cy\n", "\\end{align}\n", "\n", "Use solve to find the fixed points" ] }, { "cell_type": "code", "collapsed": false, "input": [], "language": "python", "metadata": {}, "outputs": [] } ], "metadata": {} } ] }