{ "metadata": { "name": "Lecture-5-Sympy" }, "nbformat": 3, "nbformat_minor": 0, "worksheets": [ { "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Sympy - Symbolic algebra in Python\n", "\n", "J.R. Johansson (robert@riken.jp) http://dml.riken.jp/~rob/\n", "\n", "The latest version of this [IPython notebook](http://ipython.org/notebook.html) lecture is available at [http://github.com/jrjohansson/scientific-python-lectures](http://github.com/jrjohansson/scientific-python-lectures).\n", "\n", "The other notebooks in this lecture series are indexed at [http://jrjohansson.github.com](http://jrjohansson.github.com)." ] }, { "cell_type": "code", "collapsed": false, "input": [ "%pylab inline" ], "language": "python", "metadata": {}, "outputs": [ { "output_type": "stream", "stream": "stdout", "text": [ "\n", "Welcome to pylab, a matplotlib-based Python environment [backend: module://IPython.kernel.zmq.pylab.backend_inline].\n", "For more information, type 'help(pylab)'.\n" ] } ], "prompt_number": 1 }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Introduction\n", "\n", "There are two notable Computer Algebra Systems (CAS) for Python:\n", "\n", "* [SymPy](http://sympy.org/en/index.html) - A python module that can be used in any Python program, or in an IPython session, that provides powerful CAS features. \n", "* [Sage](http://www.sagemath.org/) - Sage is a full-featured and very powerful CAS enviroment that aims to provide an open source system that competes with Mathematica and Maple. Sage is not a regular Python module, but rather a CAS environment that uses Python as its programming language.\n", "\n", "Sage is in some aspects more powerful than SymPy, but both offer very comprehensive CAS functionality. The advantage of SymPy is that it is a regular Python module and integrates well with the IPython notebook. \n", "\n", "In this lecture we will therefore look at how to use SymPy with IPython notebooks. If you are interested in an open source CAS environment I also recommend to read more about Sage.\n", "\n", "To get started using SymPy in a Python program or notebook, import the module `sympy`:" ] }, { "cell_type": "code", "collapsed": false, "input": [ "from sympy import *" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 2 }, { "cell_type": "markdown", "metadata": {}, "source": [ "To get nice-looking $\\LaTeX$ formatted output:" ] }, { "cell_type": "code", "collapsed": false, "input": [ "# commands starting with % are IPython commands\n", "%load_ext sympy.interactive.ipythonprinting\n", "\n", "# in older versions of sympy/ipython, use\n", "#%load_ext sympyprinting" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 3 }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Symbolic variables\n", "\n", "In SymPy we need to create symbols for the variables we want to work with. We can create a new symbol using the `Symbol` class:" ] }, { "cell_type": "code", "collapsed": false, "input": [ "x = Symbol('x')" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 4 }, { "cell_type": "code", "collapsed": false, "input": [ "(pi + x)**2" ], "language": "python", "metadata": {}, "outputs": [ { "latex": [ "$$\\left(x + \\pi\\right)^{2}$$" ], "metadata": {}, "output_type": "pyout", "prompt_number": 5, "text": [ " 2\n", "(x + \u03c0) " ] } ], "prompt_number": 5 }, { "cell_type": "code", "collapsed": false, "input": [ "# alternative way of defining symbols\n", "a, b, c = symbols(\"a, b, c\")" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 6 }, { "cell_type": "code", "collapsed": false, "input": [ "type(a)" ], "language": "python", "metadata": {}, "outputs": [ { "metadata": {}, "output_type": "pyout", "prompt_number": 7, "text": [ "sympy.core.symbol.Symbol" ] } ], "prompt_number": 7 }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can add assumptions to symbols when we create them:" ] }, { "cell_type": "code", "collapsed": false, "input": [ "x = Symbol('x', real=True)" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 8 }, { "cell_type": "code", "collapsed": false, "input": [ "x.is_imaginary" ], "language": "python", "metadata": {}, "outputs": [ { "latex": [ "$$False$$" ], "metadata": {}, "output_type": "pyout", "prompt_number": 9, "text": [ "False" ] } ], "prompt_number": 9 }, { "cell_type": "code", "collapsed": false, "input": [ "x = Symbol('x', positive=True)" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 10 }, { "cell_type": "code", "collapsed": false, "input": [ "x > 0" ], "language": "python", "metadata": {}, "outputs": [ { "latex": [ "$$True$$" ], "metadata": {}, "output_type": "pyout", "prompt_number": 11, "text": [ "True" ] } ], "prompt_number": 11 }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Complex numbers\n", "\n", "The imaginary unit is denoted `I` in Sympy. " ] }, { "cell_type": "code", "collapsed": false, "input": [ "1+1*I" ], "language": "python", "metadata": {}, "outputs": [ { "latex": [ "$$1 + i$$" ], "metadata": {}, "output_type": "pyout", "prompt_number": 12, "text": [ "1 + \u2148" ] } ], "prompt_number": 12 }, { "cell_type": "code", "collapsed": false, "input": [ "I**2" ], "language": "python", "metadata": {}, "outputs": [ { "latex": [ "$$-1$$" ], "metadata": {}, "output_type": "pyout", "prompt_number": 13, "text": [ "-1" ] } ], "prompt_number": 13 }, { "cell_type": "code", "collapsed": false, "input": [ "(x * I + 1)**2" ], "language": "python", "metadata": {}, "outputs": [ { "latex": [ "$$\\left(i x + 1\\right)^{2}$$" ], "metadata": {}, "output_type": "pyout", "prompt_number": 14, "text": [ " 2\n", "(\u2148\u22c5x + 1) " ] } ], "prompt_number": 14 }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Rational numbers\n", "\n", "There are three different numerical types in SymPy: `Real`, `Rational`, `Integer`: " ] }, { "cell_type": "code", "collapsed": false, "input": [ "r1 = Rational(4,5)\n", "r2 = Rational(5,4)" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 15 }, { "cell_type": "code", "collapsed": false, "input": [ "r1" ], "language": "python", "metadata": {}, "outputs": [ { "latex": [ "$$\\frac{4}{5}$$" ], "metadata": {}, "output_type": "pyout", "prompt_number": 16, "text": [ "4/5" ] } ], "prompt_number": 16 }, { "cell_type": "code", "collapsed": false, "input": [ "r1+r2" ], "language": "python", "metadata": {}, "outputs": [ { "latex": [ "$$\\frac{41}{20}$$" ], "metadata": {}, "output_type": "pyout", "prompt_number": 17, "text": [ "41\n", "\u2500\u2500\n", "20" ] } ], "prompt_number": 17 }, { "cell_type": "code", "collapsed": false, "input": [ "r1/r2" ], "language": "python", "metadata": {}, "outputs": [ { "latex": [ "$$\\frac{16}{25}$$" ], "metadata": {}, "output_type": "pyout", "prompt_number": 18, "text": [ "16\n", "\u2500\u2500\n", "25" ] } ], "prompt_number": 18 }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Numerical evaluation\n", "\n", "SymPy uses a library for artitrary precision as numerical backend, and has predefined SymPy expressions for a number of mathematical constants, such as: `pi`, `e`, `oo` for infinity.\n", "\n", "To evaluate an expression numerically we can use the `evalf` function (or `N`). It takes an argument `n` which specifies the number of significant digits." ] }, { "cell_type": "code", "collapsed": false, "input": [ "pi.evalf(n=50)" ], "language": "python", "metadata": {}, "outputs": [ { "latex": [ "$$3.1415926535897932384626433832795028841971693993751$$" ], "metadata": {}, "output_type": "pyout", "prompt_number": 19, "text": [ "3.1415926535897932384626433832795028841971693993751" ] } ], "prompt_number": 19 }, { "cell_type": "code", "collapsed": false, "input": [ "y = (x + pi)**2" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 20 }, { "cell_type": "code", "collapsed": false, "input": [ "N(y, 5) # same as evalf" ], "language": "python", "metadata": {}, "outputs": [ { "latex": [ "$$\\left(x + 3.1416\\right)^{2}$$" ], "metadata": {}, "output_type": "pyout", "prompt_number": 21, "text": [ " 2\n", "(x + 3.1416) " ] } ], "prompt_number": 21 }, { "cell_type": "markdown", "metadata": {}, "source": [ "When we numerically evaluate algebraic expressions we often want to substitute a symbol with a numerical value. In SymPy we do that using the `subs` function:" ] }, { "cell_type": "code", "collapsed": false, "input": [ "y.subs(x, 1.5)" ], "language": "python", "metadata": {}, "outputs": [ { "latex": [ "$$\\left(1.5 + \\pi\\right)^{2}$$" ], "metadata": {}, "output_type": "pyout", "prompt_number": 22, "text": [ " 2\n", "(1.5 + \u03c0) " ] } ], "prompt_number": 22 }, { "cell_type": "code", "collapsed": false, "input": [ "N(y.subs(x, 1.5))" ], "language": "python", "metadata": {}, "outputs": [ { "latex": [ "$$21.5443823618587$$" ], "metadata": {}, "output_type": "pyout", "prompt_number": 23, "text": [ "21.5443823618587" ] } ], "prompt_number": 23 }, { "cell_type": "markdown", "metadata": {}, "source": [ "The `subs` function can of course also be used to substitute Symbols and expressions:" ] }, { "cell_type": "code", "collapsed": false, "input": [ "y.subs(x, a+pi)" ], "language": "python", "metadata": {}, "outputs": [ { "latex": [ "$$\\left(a + 2 \\pi\\right)^{2}$$" ], "metadata": {}, "output_type": "pyout", "prompt_number": 24, "text": [ " 2\n", "(a + 2\u22c5\u03c0) " ] } ], "prompt_number": 24 }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can also combine numerical evolution of expressions with NumPy arrays:" ] }, { "cell_type": "code", "collapsed": false, "input": [ "import numpy" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 25 }, { "cell_type": "code", "collapsed": false, "input": [ "x_vec = numpy.arange(0, 10, 0.1)" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 26 }, { "cell_type": "code", "collapsed": false, "input": [ "y_vec = numpy.array([N(((x + pi)**2).subs(x, xx)) for xx in x_vec])" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 27 }, { "cell_type": "code", "collapsed": false, "input": [ "fig, ax = subplots()\n", "ax.plot(x_vec, y_vec);" ], "language": "python", "metadata": {}, "outputs": [ { "metadata": {}, "output_type": "display_data", "png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAEACAYAAABI5zaHAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3X9YVGXex/H3mJRtVprmoKAPrUg4iqkllRs1ptBqaWYu\nK7YrC1o9a7b92rLaH2G7CWa7qfVYa6XZjxWtLTHXyEyHNDU027SwIIMEFDYl3EUTVM7zx12kpeQM\nA2dm+Lyua64LjzPnfJtr/ezN99znvh2WZVmIiEhIaWN3ASIi4n8KdxGREKRwFxEJQQp3EZEQpHAX\nEQlBCncRkRDUaLinp6fjdDqJi4trOJafn098fDwDBgxg0KBBbNq0qeHvMjMz6dWrF7GxsaxcubL5\nqhYRkUY5GpvnvnbtWtq3b8+ECRPYtm0bAG63m/vuu4+rrrqK119/nYcffpg1a9ZQUFDA+PHj2bRp\nE+Xl5QwbNozCwkLatNEvByIiLa3R5E1ISKBjx47HHOvatSv79u0DoLq6moiICABycnJISUkhLCyM\nqKgooqOjyc/Pb6ayRUSkMW29/UBWVhaXXXYZv/3tb6mvr2fDhg0A7Nq1i0suuaThfZGRkZSXl/uv\nUhEROWle90wmTpzInDlz2LlzJ48++ijp6eknfK/D4WhScSIi4huvR+75+fmsWrUKgLFjxzJp0iQA\nIiIiKC0tbXhfWVlZQ8vmaNHR0ezYscPXekVEWqWePXvy6aefnvT7vR65R0dHk5eXB8Dq1auJiYkB\nYNSoUWRnZ1NXV0dxcTFFRUXEx8d/7/M7duzAsiy9LIsHHnjA9hoC5aXvQt+FvovGX94Oihsduaek\npJCXl8eePXvo3r07Dz74IPPmzeOWW26htraW008/nXnz5gHgcrlITk7G5XLRtm1b5s6dq7aMiIhN\nGg33RYsWHff4u+++e9zj999/P/fff3/TqxIRkSbRJHQbud1uu0sIGPouvqXv4lv6LnzX6ENMzXJB\nh4MWvqSISNDzNjs1chcRCUEKdxGREKRwFxEJQQp3EZEQpHAXEQlBCncRkQD0xBPw9WIAPlG4i4gE\nmEcegZkzoUcP38/h9cJhIiLSPCwL/vQnePFFePttiIz0/VwKdxGRAGBZcP/9sHy5CXans2nnU7iL\niNisvh7uuAPWrYM1a6Bz56afU+EuImKjI0fg5pth+3ZYvRrOPts/51W4i4jY5NAhSE2Fykp44w1o\n395/51a4i4jY4OBBGDfOBPzy5XD66f49v6ZCioi0sP37YeRIOPVUePVV/wc7/EC4p6en43Q6iYuL\nO+b4Y489Ru/evenbty9Tp05tOJ6ZmUmvXr2IjY1l5cqV/q9WRCTIVVdDUhJ07w6LFpmAbw6NtmXS\n0tK49dZbmTBhQsOxNWvWsGzZMrZu3UpYWBhffPEFAAUFBSxevJiCggLKy8sZNmwYhYWFtGmjXw5E\nRAD+/W/46U/hsstg1ixoznhs9NQJCQl07NjxmGNPPPEE9913H2FhYQCce+65AOTk5JCSkkJYWBhR\nUVFER0eTn5/fTGWLiASX0lK4/HK4+mqYPbt5gx186LkXFRXx9ttvc8kll+B2u9m8eTMAu3btIvKo\nx6kiIyMpLy/3X6UiIkGqqAgSEmDSJPMEqsPR/Nf0erbM4cOH+fLLL9m4cSObNm0iOTmZzz777Ljv\ndZzgvyAjI6PhZ7fbrX0SRSRkbd0Kw4dDRgbceOPJf87j8eDxeHy+rtfhHhkZyZgxYwAYNGgQbdq0\nYc+ePURERFBaWtrwvrKyMiIiIo57jqPDXUQkVK1fD9ddB3PmwM9/7t1nvzvwnTZtmlef97otM3r0\naFavXg1AYWEhdXV1dO7cmVGjRpGdnU1dXR3FxcUUFRURHx/v7elFRELCypVw7bWwcKH3we4PjY7c\nU1JSyMvLY+/evXTv3p0HH3yQ9PR00tPTiYuL49RTT+W5554DwOVykZycjMvlom3btsydO/eEbRkR\nkVD28sswebKZw37ZZfbU4LAsy2rRCzoctPAlRURazFNPwQMPwIoV0L+//87rbXZq+QERET+ZMQOe\nfNLsoNSrl721KNxFRJrIsmDqVDNaX7cOTjCXpEUp3EVEmuDwYbjpJrNkb14edOpkd0WGwl1ExEdf\nfWVWdqythVWr4Iwz7K7oW1r4RUTEB/v2mXVizjgDli0LrGAHhbuIiNd274YrroB+/eCFF5pvZcem\nULiLiHjh00/N3PWxY82Tp4G68K167iIiJ2nLFrjmGpg2zbt1YuygcBcROQmrVsH48fC3v5n1YgJd\ngP5CISISOLKzTbC//HJwBDto5C4i0qjZs2HmTHjrLfjOjqMBTeEuInIc9fVw771mmuO6dRAVZXdF\n3lG4i4h8R10dpKfDjh3wzjuB89SpNxTuIiJH+e9/4frroV0704r50Y/srsg3uqEqIvK1bx5OioqC\nV14J3mAHhbuICAAffwyDB8OYMWa6Y9sg72s0Gu7p6ek4nU7ijnOL+C9/+Qtt2rShqqqq4VhmZia9\nevUiNjaWlStX+r9aEZFmsH49uN3wxz/C738PobCJXKPhnpaWRm5u7veOl5aW8uabb/I///M/DccK\nCgpYvHgxBQUF5ObmMnnyZOrr6/1fsYiIH73yitnrdMECSEuzuxr/aTTcExIS6Nix4/eO33nnnTz8\n8MPHHMvJySElJYWwsDCioqKIjo4mPz/fv9WKiPjR7Nlw663wxhswfLjd1fiX112lnJwcIiMj6dev\n3zHHd+3axSWXXNLw58jISMrLy5teoYiIn9XXw913m52T3nkn+Oawnwyvwv3AgQNMnz6dN998s+FY\nYxu2Ok7QuMrIyGj42e1243a7vSlDRMRnX30FEyZAZaUJ9nPOsbui4/N4PHg8Hp8/71W479ixg5KS\nEi644AIAysrKuPDCC3n33XeJiIigtLS04b1lZWVEnGAjwaPDXUSkpezZY/rrPXrAypVmLnug+u7A\nd9q0aV593qupkHFxcVRWVlJcXExxcTGRkZFs2bIFp9PJqFGjyM7Opq6ujuLiYoqKioiPj/eqGBGR\n5vLpp2aqY0ICvPhiYAe7PzQa7ikpKQwePJjCwkK6d+/OggULjvn7o9suLpeL5ORkXC4Xw4cPZ+7c\nuSdsy4iItKT1680GG3feCVlZgbvBhj85rMaa5s1xQYej0T69iIg/vfQSTJ4Mzz0X3DNivM3OIH8G\nS0Tk+CzLLNX72GPw5pvQv7/dFbUshbuIhJxDh+CWWyA/HzZsgMhIuytqeQp3EQkp+/bBz34GYWGw\ndi2ceabdFdmjFdxWEJHWoqQEfvITiImBnJzWG+ygcBeRELFxo5nqeNNNps8e7Ks6NlUr/88XkVCw\nZInpsS9YANdcY3c1gUHhLiJBy7LgoYdg3jxYtQq+fnheULiLSJCqrYVJk+CTT+Ddd6FrV7srCizq\nuYtI0PniCxg6FA4eBI9HwX48CncRCSoffgjx8WbnpMWLg3uf0+aktoyIBI1//tPslvToo3DDDXZX\nE9gU7iIS8CzLBPojj5j565deandFgU/hLiIBrbYWfv1reO89s5TAUVs3SyMU7iISsL74AsaMgc6d\nza5J7dvbXVHw0A1VEQlI27aZG6dXXAH/+IeC3VsauYtIwFm61CwjMGsWjB9vdzXBqdGRe3p6Ok6n\nk7i4uIZjd999N7179+aCCy5gzJgx7Nu3r+HvMjMz6dWrF7GxsaxcubL5qhaRkGRZ8Oc/w623mpkx\nCnbfNRruaWlp5ObmHnMsKSmJjz76iA8++ICYmBgyMzMBKCgoYPHixRQUFJCbm8vkyZOpr69vvspF\nJKTs3w8//zksX27WYR80yO6Kgluj4Z6QkEDHjh2POZaYmEibrzcgvPjiiykrKwMgJyeHlJQUwsLC\niIqKIjo6mvz8/GYqW0RCyeefm6V6Tz9dT5z6S5NuqM6fP58RI0YAsGvXLiKP2u4kMjKS8vLyplUn\nIiEvLw8uuQRSU+HZZ6FdO7srCg0+31B96KGHOPXUUxnfSFPM4XAc93hGRkbDz263G7fb7WsZIhKk\nLAvmzoUHH4Tnn4ekJLsrCiwejwePx+Pz530K92effZYVK1bw1ltvNRyLiIigtLS04c9lZWVEREQc\n9/NHh7uItD61tWb99Y0bYf166NnT7ooCz3cHvtOmTfPq8163ZXJzc5k5cyY5OTm0O+r3p1GjRpGd\nnU1dXR3FxcUUFRURHx/v7elFJMTt3g1DhkBVlXniVMHePBoN95SUFAYPHswnn3xC9+7dmT9/Prfe\neis1NTUkJiYyYMAAJk+eDIDL5SI5ORmXy8Xw4cOZO3fuCdsyItI6rV9vZsEMHw4vv9y69zhtbg7L\nsqwWvaDDQQtfUkQCwLx58Pvfw/z52grPF95mp55QFZFmVVtrHkpat868YmLsrqh10NoyItJsysrM\n2jB795qt8BTsLUfhLiLN4u23zcJfo0erv24HtWVExK8sC2bPhsxMzV+3k8JdRPympgZuvBE++cS0\nYaKi7K6o9VJbRkT8oqjILCPQrp3ZWEPBbi+Fu4g02auvmoW/pkwxUx1PP93uikRtGRHx2eHD8Lvf\nQXa2WapXD6UHDoW7iPikshJSUuCUU8zm1Z07212RHE1tGRHx2rp1cOGFphWTm6tgD0QauYvISbMs\nePRRmDHDrL0+fLjdFcmJKNxF5KRUV0N6OpSWappjMFBbRkR+0Pvvw0UXQbdupiWjYA98CncROSHL\nMqs5JiXBQw/B44/DaafZXZWcDLVlROS4amrg5pth61ZYuxZiY+2uSLzR6Mg9PT0dp9NJXFxcw7Gq\nqioSExOJiYkhKSmJ6urqhr/LzMykV69exMbGsnLlyuarWkSa1bZtpg3Trp3pryvYg0+j4Z6WlkZu\nbu4xx7KyskhMTKSwsJChQ4eSlZUFQEFBAYsXL6agoIDc3FwmT55MfX1981UuIn5nWfDMM3DllXDf\nfebnH/3I7qrEF42Ge0JCAh07djzm2LJly0hNTQUgNTWVpUuXApCTk0NKSgphYWFERUURHR1Nfn5+\nM5UtIv5WUwO//KWZ6piXB1//M5cg5fUN1crKSpxOJwBOp5PKykoAdu3aRWRkZMP7IiMjKS8v91OZ\nItKcPvjAtGFOOw3y88HlsrsiaaomzZZxOByNboKtDbJFAptlwRNPwNChZo0YtWFCh9ezZZxOJxUV\nFYSHh7N79266dOkCQEREBKWlpQ3vKysrIyIi4rjnyMjIaPjZ7Xbjdru9LUNEmqi62qy9XlRklug9\n/3y7K5KjeTwePB6Pz593WD+wnXZJSQkjR45k27ZtANxzzz106tSJqVOnkpWVRXV1NVlZWRQUFDB+\n/Hjy8/MpLy9n2LBhfPrpp98bvXu7g7eI+N/GjWbRrxEj4C9/MbNiJLB5m52NjtxTUlLIy8tjz549\ndO/enQcffJB7772X5ORknnnmGaKioliyZAkALpeL5ORkXC4Xbdu2Ze7cuWrLiASY+nqYORP++ld4\n8km47jq7K5Lm8oMjd79fUCN3EVtUVMCECXDgAPz979Cjh90ViTe8zU4tPyDSCqxYAQMGwKWXgsej\nYG8NtPyASAirrTUPI738stkt6Yor7K5IWorCXSREbd9ubpqed55Z1bFTJ7srkpaktoxIiPlmJceE\nBJg8GV55RcHeGmnkLhJC9uwxc9eLi81Kjr17212R2EUjd5EQ8eab0L8/9OxpVnJUsLduGrmLBLmD\nB83SAYsXm31Nhw2zuyIJBAp3kSC2bRvccAP06mUW/1JvXb6htoxIEKqvN0vzXnkl3HmnmeqoYJej\naeQuEmRKS+FXv4KvvjK99R//2O6KJBBp5C4SRP7+d7jwQrNE79tvK9jlxDRyFwkCe/fCLbeYzapz\nc2HgQLsrkkCnkbtIgHv9dbjgAujaFd57T8EuJ0cjd5EAVVMDd99twv3552HIELsrkmCikbtIAFq7\n1ozWDx40UxwV7OItjdxFAshXX8Ef/mBunD75JIwaZXdFEqx8HrlnZmbSp08f4uLiGD9+PLW1tVRV\nVZGYmEhMTAxJSUlUV1f7s1aRkPbuu6afvnOnuXGqYJem8CncS0pKeOqpp9iyZQvbtm3jyJEjZGdn\nk5WVRWJiIoWFhQwdOpSsrCx/1ysScmpr4f774dprYdo0WLIEOne2uyoJdj6F+1lnnUVYWBgHDhzg\n8OHDHDhwgG7durFs2TJSU1MBSE1NZenSpX4tViTUbN5s5q1v325668nJdlckocKncD/nnHO46667\n6NGjB926daNDhw4kJiZSWVmJ0+kEwOl0UllZ6ddiRUJFba1Z7Ovqq82o/ZVX4Ot/OiJ+4VO479ix\ng1mzZlFSUsKuXbuoqanhhRdeOOY9DocDh8PhlyJFQkl+vhmtf/SRGa2PHw/6pyL+5tNsmc2bNzN4\n8GA6fb1S0ZgxY9iwYQPh4eFUVFQQHh7O7t276dKly3E/n5GR0fCz2+3G7Xb7UoZIUPnqK3jgAXju\nObPo17hxCnU5MY/Hg8fj8fnzDsuyLG8/9MEHH3DDDTewadMm2rVrx69+9Svi4+P5/PPP6dSpE1On\nTiUrK4vq6urv3VR1OBz4cEmRoLZuHUycaDbTeOwxOMG4R+SEvM1On8Id4OGHH2bhwoW0adOGgQMH\n8vTTT/Pf//6X5ORkdu7cSVRUFEuWLKFDhw5NKlAkmNXUwH33mZ76Y4/BmDF2VyTBqsXC3VcKd2kt\nVq6Em28Gtxv++lfo2NHuiiSYeZudekJVxM/27jUbaOTlwd/+BlddZXdF0hppbRkRP7EsyM6Gvn2h\nQwf48EMFu9hHI3cRP9i5EyZPhpISePVVuOQSuyuS1k4jd5EmOHLE3CgdONAE+pYtCnYJDBq5i/jo\ngw/gppugXTt45x04/3y7KxL5lkbuIl7avx/uuQcSE024r1mjYJfAo3AX8cKKFeaGaXm5uWE6cSK0\n0b8iCUBqy4ichF274Lbb4P33zfTGpCS7KxJpnMYcIo04fBjmzIF+/SA2FrZtU7BLcNDIXeQE8vPh\nf/8Xzj7b7Gnau7fdFYmcPI3cRb6jqgp+/WuzM9Kdd8Lq1Qp2CT4Kd5Gv1dfDggXgcsEpp5jdkX7x\nCy3LK8FJbRkRzJz1W26Bujr45z/NZhoiwUwjd2nVqqvhN78xN0knTIANGxTsEhoU7tIqfdOC6d3b\n7GdaUGAeSDrlFLsrE/EPtWWk1dm8GaZMMT+/9hpcdJG99Yg0B59H7tXV1YwdO5bevXvjcrl49913\nqaqqIjExkZiYGJKSkqiurvZnrSJN8u9/w403wsiRZhON9esV7BK6fA732267jREjRrB9+3a2bt1K\nbGwsWVlZJCYmUlhYyNChQ7+3f6qIHQ4dgtmzoU8fOPNMMwsmLU3LBkho82mbvX379jFgwAA+++yz\nY47HxsaSl5eH0+mkoqICt9vNxx9/fOwFtc2etKCVK+H22yEyEmbNMtMcRYKRt9np09iluLiYc889\nl7S0NAYOHMiNN97I/v37qaysxOl0AuB0OqmsrPTl9CJNVlQEo0aZDTSysuCNNxTs0rr4dEP18OHD\nbNmyhccff5xBgwZx++23f68F43A4cJzg6Y+MjIyGn91uN26325cyRL6nuhr+/Gd49lmYOhVeeglO\nO83uqkS85/F48Hg8Pn/ep7ZMRUUFl156KcXFxQCsW7eOzMxMPvvsM9asWUN4eDi7d+9myJAhastI\nizh8GJ5+GjIyzIj9T3+Cr3+JFAkJLdKWCQ8Pp3v37hQWFgKwatUq+vTpw8iRI1m4cCEACxcuZPTo\n0b6cXsQrublwwQWweLH5ed48BbuITyN3gA8++IBJkyZRV1dHz549WbBgAUeOHCE5OZmdO3cSFRXF\nkiVL6NChw7EX1Mhd/OTDD+G3v4XPPoNHHjFTHLUOjIQqb7PT53D3lcJdmqqiAv74R1i6FH73O7OC\n46mn2l2VSPNqkbaMiB327ze99L59zRrrn3xidkdSsIt8n8JdAt43N0tjYuCjj2DTJpg5Ezp2tLsy\nkcCltWUkYFmWWX536lQ491zThhk0yO6qRIKDwl0C0vr1JtSrqsxDSNdco5ulIt5QW0YCyvbtcN11\nMG4cpKfD1q2aBSPiC4W7BISdO02YX345DB5sbpampWl9dRFfKdzFVl98YTah7t8fwsPNmjB33w2n\nn253ZSLBTeEutti3z8xVj401+5Z+9BFMnw7feeZNRHykcJcWtX8/zJgBvXrB55+bXZEefxy6drW7\nMpHQotky0iIOHoSnnoLMTLjsMsjLM/uXikjzULhLs6qrg/nz4aGHzOJe//wnDBhgd1UioU/hLs3i\n0CF47jmztvr558M//gHx8XZXJdJ6KNzFrw4dguefN6Hesye88AL85Cd2VyXS+ijcxS++Gak/9BBE\nRcHChZCQYHdVIq2Xwl2apK7OBPn06WakrlAXCQwKd/HJwYPwzDNmWmPv3mq/iASaJs1zP3LkCAMG\nDGDkyJEAVFVVkZiYSExMDElJSVRXV/ulSAkc+/fDX/9qRum5ufDyy/DGGwp2kUDTpHCfPXs2LpcL\nx9erOmVlZZGYmEhhYSFDhw4lKyvLL0WK/aqrzU3S886DDRtg+XJ47TXNgBEJVD6He1lZGStWrGDS\npEkNWz8tW7aM1NRUAFJTU1m6dKl/qhTbVFbCvfeakXpREbz9Nrz0kuaqiwQ6n8P9jjvuYObMmbRp\n8+0pKisrcX697bzT6aSysrLpFYotPvsMbrnF9NNrauC998zN0thYuysTkZPh0w3V5cuX06VLFwYM\nGIDH4znuexwOR0O75rsyMjIafna73bjdbl/KkGbwr3+Zm6Rvvgk33WTWV//6/69FpAV5PJ4T5uvJ\ncFjebKf9tfvvv5/nn3+etm3bcvDgQf7zn/8wZswYNm3ahMfjITw8nN27dzNkyBA+/vjjYy/o5Q7e\n0vwsC1avhocfhg8/hDvuMMF+1ll2VyYi3/A2O30K96Pl5eXxyCOP8Nprr3HPPffQqVMnpk6dSlZW\nFtXV1d+7qapwDxyHDpnZLjNnmqmNd98N48fDaafZXZmIfJe32emXee7ftF/uvfdekpOTeeaZZ4iK\nimLJkiX+OL342X/+A08/DbNmmdkv06bB1VdDGy0ALRIymjxy9/qCGrnb5vPPYc4cePZZSEyEu+6C\nQYPsrkpEToa32amxWoizLDMv/ec/h4EDzUbTW7ZAdraCXSSUafmBEPVNP33WLLNP6W9+YzbL0E1S\nkdZBbZkQ8+9/w7x58OSTZiu722+Ha66BU06xuzIRaQq1ZVqpzZshNdVsjPH557BiBaxZA9deq2AX\naY00cg9itbVmKYDHH4eKCvj1r2HSJOjUye7KRMTfWnyeu7cU7k1XUgJ/+5vZm7RfP5gyRa0XkVCn\ntkyIOnLEbC49ciRceKF56Ojtt80yAWq9iMh3abZMgNu924zQ580za7zcfDMsXgw/+pHdlYlIIFO4\nB6D6eli50gT6mjXws5/BK6+YEbuIyMlQuAeQsjIzSp8/39wUvflms8zumWfaXZmIBBuFu83q6syO\nRvPnw8aNMG4cvPqqNsMQkaZRuNtk2zZYsMBsLO1yQXq6mdaoXrqI+IPCvQXt3QuLFplQr6yECRNg\n/XqIjra7MhEJNZrn3swOHTJPiy5cCG+9BSNGQFoaDB2q6YsicvL0EFMAsCzIzzctl8WLzZIAqalm\n1svZZ9tdnYgEI1s26xDj00/hxRfNC+CXvzQ3SX/8Y3vrEpHWx6cnVEtLSxkyZAh9+vShb9++zJkz\nB4CqqioSExOJiYkhKSmJ6upqvxYbiHbvhtmz4eKL4Sc/MX3155+HTz6BP/xBwS4i9vCpLVNRUUFF\nRQX9+/enpqaGCy+8kKVLl7JgwQI6d+7MPffcw4wZM/jyyy9Dcg/VqiozXXHRInjvPbMkwPjxMGwY\ntNXvQiLSDGzpuY8ePZopU6YwZcoU8vLycDqdVFRU4Ha7+fjjj5tUYKCoroZly0wPfd06s01dSoq5\nQXr66XZXJyKhrsXDvaSkhCuuuIIPP/yQHj168OWXXwJgWRbnnHNOw599LdBO3wT6Sy9BXh4MGWK2\nqxs5Uk+NikjLatEbqjU1NVx//fXMnj2bM7+Tdg6HA4fDcdzPZWRkNPzsdrtxu91NKcOv9uyBnByz\nRd0778CVV5qnRl98UVvUiUjL8Xg8eDwenz/v88j90KFDXHPNNQwfPpzbb78dgNjYWDweD+Hh4eze\nvZshQ4YERVumtBSWLjV99Pfeg6QkGDvWtFw0QheRQNAi67lblsXEiRNxuVwNwQ4watQoFi5cCMDC\nhQsZPXq0L6dvdpZlHv//859h0CDo39+E+m23mR2NXnrJtF8U7CISrHwaua9bt47LL7+cfv36NbRe\nMjMziY+PJzk5mZ07dxIVFcWSJUvo0KHDsRe0aeR+6BCsXWt66MuWmWV1r73WvBISICysxUsSETlp\nekL1KHv2wOuvw/LlZsei6GgYNcrcEO3XD05wS0BEJOC06nCvr4f33zdruaxYAQUFZg2Xq682/fOu\nXZvlsiIiza7VhfsXX5hReW4uvPEGdOxognzECNNuOe00v11KRMQ2IR/utbWwYYPZhm7lSigqMvPP\nf/pTuOoqOO88PxYrIhIgQi7c6+th61azXO6qVWbu+fnnmyBPSoJLL9XNUBEJfUEf7pZlFt3yeEyg\nr1kD55xj1m0ZNgzcbvNnEZHWJOjC3bLg44/N4/15eSbUTz3VhPjQoabl0r17S1YoIhJ4giLcN22y\nWLvWLMC1di20bw+XXw5XXGHCPCqqJSsSEQl8QRHuffpYJCTAZZeZUNfIXESkcUER7oG2toyISKBr\nkbVlREQksCncRURCkMJdRCQEKdxFREKQwl1EJAQp3EVEQpDfwz03N5fY2Fh69erFjBkz/H16ERE5\nCX4N9yNHjjBlyhRyc3MpKChg0aJFbN++3Z+XCClN2fw21Oi7+Ja+i2/pu/CdX8M9Pz+f6OhooqKi\nCAsLY9y4ceTk5PjzEiFF/8P9lr6Lb+m7+Ja+C9/5NdzLy8vpftRaApGRkZSXl/vzEiIichL8Gu4O\nbUoqIhKh55+5AAAEZ0lEQVQYLD/asGGDddVVVzX8efr06VZWVtYx7+nZs6cF6KWXXnrp5cWrZ8+e\nXuWxXxcOO3z4MOeffz5vvfUW3bp1Iz4+nkWLFtG7d29/XUJERE5CW7+erG1bHn/8ca666iqOHDnC\nxIkTFewiIjZo8SV/RUSk+bXoE6p6wMkoLS1lyJAh9OnTh759+zJnzhy7S7LdkSNHGDBgACNHjrS7\nFFtVV1czduxYevfujcvlYuPGjXaXZJvMzEz69OlDXFwc48ePp7a21u6SWkx6ejpOp5O4uLiGY1VV\nVSQmJhITE0NSUhLV1dWNnqPFwl0POH0rLCyMRx99lI8++oiNGzfyf//3f632u/jG7NmzcblcrX7G\n1W233caIESPYvn07W7dubbVtzZKSEp566im2bNnCtm3bOHLkCNnZ2XaX1WLS0tLIzc095lhWVhaJ\niYkUFhYydOhQsrKyGj1Hi4W7HnD6Vnh4OP379wegffv29O7dm127dtlclX3KyspYsWIFkyZNatW7\ndO3bt4+1a9eSnp4OmHtYZ599ts1V2eOss84iLCyMAwcOcPjwYQ4cOEBERITdZbWYhIQEOnbseMyx\nZcuWkZqaCkBqaipLly5t9BwtFu56wOn4SkpKeP/997n44ovtLsU2d9xxBzNnzqRNm9a9jl1xcTHn\nnnsuaWlpDBw4kBtvvJEDBw7YXZYtzjnnHO666y569OhBt27d6NChA8OGDbO7LFtVVlbidDoBcDqd\nVFZWNvr+FvvX1Np/3T6empoaxo4dy+zZs2nfvr3d5dhi+fLldOnShQEDBrTqUTuYqcRbtmxh8uTJ\nbNmyhTPOOOMHf/UOVTt27GDWrFmUlJSwa9cuampqePHFF+0uK2A4HI4fzNQWC/eIiAhKS0sb/lxa\nWkpkZGRLXT7gHDp0iOuvv55f/OIXjB492u5ybLN+/XqWLVvGeeedR0pKCqtXr2bChAl2l2WLyMhI\nIiMjGTRoEABjx45ly5YtNldlj82bNzN48GA6depE27ZtGTNmDOvXr7e7LFs5nU4qKioA2L17N126\ndGn0/S0W7hdddBFFRUWUlJRQV1fH4sWLGTVqVEtdPqBYlsXEiRNxuVzcfvvtdpdjq+nTp1NaWkpx\ncTHZ2dlceeWVPPfcc3aXZYvw8HC6d+9OYWEhAKtWraJPnz42V2WP2NhYNm7cyFdffYVlWaxatQqX\ny2V3WbYaNWoUCxcuBGDhwoU/PChs0noDXlqxYoUVExNj9ezZ05o+fXpLXjqgrF271nI4HNYFF1xg\n9e/f3+rfv7/1+uuv212W7TwejzVy5Ei7y7DVv/71L+uiiy6y+vXrZ1133XVWdXW13SXZZsaMGZbL\n5bL69u1rTZgwwaqrq7O7pBYzbtw4q2vXrlZYWJgVGRlpzZ8/39q7d681dOhQq1evXlZiYqL15Zdf\nNnoOPcQkIhKCWvf0BBGREKVwFxEJQQp3EZEQpHAXEQlBCncRkRCkcBcRCUEKdxGREKRwFxEJQf8P\naFlJKJFyhJ8AAAAASUVORK5CYII=\n", "text": [ "" ] } ], "prompt_number": 28 }, { "cell_type": "markdown", "metadata": {}, "source": [ "However, this kind of numerical evolution can be very slow, and there is a much more efficient way to do it: Use the function `lambdify` to \"compile\" a Sympy expression into a function that is much more efficient to evaluate numerically:" ] }, { "cell_type": "code", "collapsed": false, "input": [ "f = lambdify([x], (x + pi)**2, 'numpy') # the first argument is a list of variables that f will be a function of: in this case only x -> f(x)" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 29 }, { "cell_type": "code", "collapsed": false, "input": [ "y_vec = f(x_vec) # now we can directly pass a numpy array and f(x) is efficiently evaluated" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 30 }, { "cell_type": "markdown", "metadata": {}, "source": [ "The speedup when using \"lambdified\" functions instead of direct numerical evaluation can be significant, often several orders of magnitude. Even in this simple example we get a significant speed up:" ] }, { "cell_type": "code", "collapsed": false, "input": [ "%%timeit\n", "\n", "y_vec = numpy.array([N(((x + pi)**2).subs(x, xx)) for xx in x_vec])" ], "language": "python", "metadata": {}, "outputs": [ { "output_type": "stream", "stream": "stdout", "text": [ "10 loops, best of 3: 20.4 ms per loop\n" ] } ], "prompt_number": 31 }, { "cell_type": "code", "collapsed": false, "input": [ "%%timeit\n", "\n", "y_vec = f(x_vec)" ], "language": "python", "metadata": {}, "outputs": [ { "output_type": "stream", "stream": "stdout", "text": [ "100000 loops, best of 3: 3.67 \u00b5s per loop\n" ] } ], "prompt_number": 32 }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Algebraic manipulations\n", "\n", "One of the main uses of an CAS is to perform algebraic manipulations of expressions. For example, we might want to expand a product, factor an expression, or simply an expression. The functions for doing these basic operations in SymPy are demonstrated in this section." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Expand and factor\n", "\n", "The first steps in an algebraic manipulation " ] }, { "cell_type": "code", "collapsed": false, "input": [ "(x+1)*(x+2)*(x+3)" ], "language": "python", "metadata": {}, "outputs": [ { "latex": [ "$$\\left(x + 1\\right) \\left(x + 2\\right) \\left(x + 3\\right)$$" ], "metadata": {}, "output_type": "pyout", "prompt_number": 33, "text": [ "(x + 1)\u22c5(x + 2)\u22c5(x + 3)" ] } ], "prompt_number": 33 }, { "cell_type": "code", "collapsed": false, "input": [ "expand((x+1)*(x+2)*(x+3))" ], "language": "python", "metadata": {}, "outputs": [ { "latex": [ "$$x^{3} + 6 x^{2} + 11 x + 6$$" ], "metadata": {}, "output_type": "pyout", "prompt_number": 34, "text": [ " 3 2 \n", "x + 6\u22c5x + 11\u22c5x + 6" ] } ], "prompt_number": 34 }, { "cell_type": "markdown", "metadata": {}, "source": [ "The `expand` function takes a number of keywords arguments which we can tell the functions what kind of expansions we want to have performed. For example, to expand trigonometric expressions, use the `trig=True` keyword argument:" ] }, { "cell_type": "code", "collapsed": false, "input": [ "sin(a+b)" ], "language": "python", "metadata": {}, "outputs": [ { "latex": [ "$$\\sin{\\left (a + b \\right )}$$" ], "metadata": {}, "output_type": "pyout", "prompt_number": 35, "text": [ "sin(a + b)" ] } ], "prompt_number": 35 }, { "cell_type": "code", "collapsed": false, "input": [ "expand(sin(a+b), trig=True)" ], "language": "python", "metadata": {}, "outputs": [ { "latex": [ "$$\\sin{\\left (a \\right )} \\cos{\\left (b \\right )} + \\sin{\\left (b \\right )} \\cos{\\left (a \\right )}$$" ], "metadata": {}, "output_type": "pyout", "prompt_number": 36, "text": [ "sin(a)\u22c5cos(b) + sin(b)\u22c5cos(a)" ] } ], "prompt_number": 36 }, { "cell_type": "markdown", "metadata": {}, "source": [ "See `help(expand)` for a detailed explanation of the various types of expansions the `expand` functions can perform." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The opposite a product expansion is of course factoring. The factor an expression in SymPy use the `factor` function: " ] }, { "cell_type": "code", "collapsed": false, "input": [ "factor(x**3 + 6 * x**2 + 11*x + 6)" ], "language": "python", "metadata": {}, "outputs": [ { "latex": [ "$$\\left(x + 1\\right) \\left(x + 2\\right) \\left(x + 3\\right)$$" ], "metadata": {}, "output_type": "pyout", "prompt_number": 37, "text": [ "(x + 1)\u22c5(x + 2)\u22c5(x + 3)" ] } ], "prompt_number": 37 }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Simplify\n", "\n", "The `simplify` tries to simplify an expression into a nice looking expression, using various techniques. More specific alternatives to the `simplify` functions also exists: `trigsimp`, `powsimp`, `logcombine`, etc. \n", "\n", "The basic usages of these functions are as follows:" ] }, { "cell_type": "code", "collapsed": false, "input": [ "# simplify expands a product\n", "simplify((x+1)*(x+2)*(x+3))" ], "language": "python", "metadata": {}, "outputs": [ { "latex": [ "$$\\left(x + 1\\right) \\left(x + 2\\right) \\left(x + 3\\right)$$" ], "metadata": {}, "output_type": "pyout", "prompt_number": 38, "text": [ "(x + 1)\u22c5(x + 2)\u22c5(x + 3)" ] } ], "prompt_number": 38 }, { "cell_type": "code", "collapsed": false, "input": [ "# simplify uses trigonometric identities\n", "simplify(sin(a)**2 + cos(a)**2)" ], "language": "python", "metadata": {}, "outputs": [ { "latex": [ "$$1$$" ], "metadata": {}, "output_type": "pyout", "prompt_number": 39, "text": [ "1" ] } ], "prompt_number": 39 }, { "cell_type": "code", "collapsed": false, "input": [ "simplify(cos(x)/sin(x))" ], "language": "python", "metadata": {}, "outputs": [ { "latex": [ "$$\\frac{1}{\\tan{\\left (x \\right )}}$$" ], "metadata": {}, "output_type": "pyout", "prompt_number": 40, "text": [ " 1 \n", "\u2500\u2500\u2500\u2500\u2500\u2500\n", "tan(x)" ] } ], "prompt_number": 40 }, { "cell_type": "markdown", "metadata": {}, "source": [ "## apart and together\n", "\n", "To manipulate symbolic expressions of fractions, we can the `apart` and `together` functions:" ] }, { "cell_type": "code", "collapsed": false, "input": [ "f1 = 1/((a+1)*(a+2))" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 41 }, { "cell_type": "code", "collapsed": false, "input": [ "f1" ], "language": "python", "metadata": {}, "outputs": [ { "latex": [ "$$\\frac{1}{\\left(a + 1\\right) \\left(a + 2\\right)}$$" ], "metadata": {}, "output_type": "pyout", "prompt_number": 42, "text": [ " 1 \n", "\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\n", "(a + 1)\u22c5(a + 2)" ] } ], "prompt_number": 42 }, { "cell_type": "code", "collapsed": false, "input": [ "apart(f1)" ], "language": "python", "metadata": {}, "outputs": [ { "latex": [ "$$- \\frac{1}{a + 2} + \\frac{1}{a + 1}$$" ], "metadata": {}, "output_type": "pyout", "prompt_number": 43, "text": [ " 1 1 \n", "- \u2500\u2500\u2500\u2500\u2500 + \u2500\u2500\u2500\u2500\u2500\n", " a + 2 a + 1" ] } ], "prompt_number": 43 }, { "cell_type": "code", "collapsed": false, "input": [ "f2 = 1/(a+2) + 1/(a+3)" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 44 }, { "cell_type": "code", "collapsed": false, "input": [ "f2" ], "language": "python", "metadata": {}, "outputs": [ { "latex": [ "$$\\frac{1}{a + 3} + \\frac{1}{a + 2}$$" ], "metadata": {}, "output_type": "pyout", "prompt_number": 45, "text": [ " 1 1 \n", "\u2500\u2500\u2500\u2500\u2500 + \u2500\u2500\u2500\u2500\u2500\n", "a + 3 a + 2" ] } ], "prompt_number": 45 }, { "cell_type": "code", "collapsed": false, "input": [ "together(f2)" ], "language": "python", "metadata": {}, "outputs": [ { "latex": [ "$$\\frac{2 a + 5}{\\left(a + 2\\right) \\left(a + 3\\right)}$$" ], "metadata": {}, "output_type": "pyout", "prompt_number": 46, "text": [ " 2\u22c5a + 5 \n", "\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\n", "(a + 2)\u22c5(a + 3)" ] } ], "prompt_number": 46 }, { "cell_type": "markdown", "metadata": {}, "source": [ "Simplify usually combines fractions but does not factor: " ] }, { "cell_type": "code", "collapsed": false, "input": [ "simplify(f2)" ], "language": "python", "metadata": {}, "outputs": [ { "latex": [ "$$\\frac{2 a + 5}{\\left(a + 2\\right) \\left(a + 3\\right)}$$" ], "metadata": {}, "output_type": "pyout", "prompt_number": 47, "text": [ " 2\u22c5a + 5 \n", "\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\n", "(a + 2)\u22c5(a + 3)" ] } ], "prompt_number": 47 }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Calculus\n", "\n", "In addition to algebraic manipulations, the other main use of CAS is to do calculus, like derivatives and integrals of algebraic expressions." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Differentiation\n", "\n", "Differentiation is usually simple. Use the `diff` function. The first argument is the expression to take the derivative of, and the second argument is the symbol by which to take the deriative:" ] }, { "cell_type": "code", "collapsed": false, "input": [ "y" ], "language": "python", "metadata": {}, "outputs": [ { "latex": [ "$$\\left(x + \\pi\\right)^{2}$$" ], "metadata": {}, "output_type": "pyout", "prompt_number": 48, "text": [ " 2\n", "(x + \u03c0) " ] } ], "prompt_number": 48 }, { "cell_type": "code", "collapsed": false, "input": [ "diff(y**2, x)" ], "language": "python", "metadata": {}, "outputs": [ { "latex": [ "$$4 \\left(x + \\pi\\right)^{3}$$" ], "metadata": {}, "output_type": "pyout", "prompt_number": 49, "text": [ " 3\n", "4\u22c5(x + \u03c0) " ] } ], "prompt_number": 49 }, { "cell_type": "markdown", "metadata": {}, "source": [ "For higher order derivatives we can do:" ] }, { "cell_type": "code", "collapsed": false, "input": [ "diff(y**2, x, x)" ], "language": "python", "metadata": {}, "outputs": [ { "latex": [ "$$12 \\left(x + \\pi\\right)^{2}$$" ], "metadata": {}, "output_type": "pyout", "prompt_number": 50, "text": [ " 2\n", "12\u22c5(x + \u03c0) " ] } ], "prompt_number": 50 }, { "cell_type": "code", "collapsed": false, "input": [ "diff(y**2, x, 2) # same as above" ], "language": "python", "metadata": {}, "outputs": [ { "latex": [ "$$12 \\left(x + \\pi\\right)^{2}$$" ], "metadata": {}, "output_type": "pyout", "prompt_number": 51, "text": [ " 2\n", "12\u22c5(x + \u03c0) " ] } ], "prompt_number": 51 }, { "cell_type": "markdown", "metadata": {}, "source": [ "To calculate the derivative of a multivariate expression, we can do:" ] }, { "cell_type": "code", "collapsed": false, "input": [ "x, y, z = symbols(\"x,y,z\")" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 52 }, { "cell_type": "code", "collapsed": false, "input": [ "f = sin(x*y) + cos(y*z)" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 53 }, { "cell_type": "markdown", "metadata": {}, "source": [ "$\\frac{d^3f}{dxdy^2}$" ] }, { "cell_type": "code", "collapsed": false, "input": [ "diff(f, x, 1, y, 2)" ], "language": "python", "metadata": {}, "outputs": [ { "latex": [ "$$- x \\left(x y \\cos{\\left (x y \\right )} + 2 \\sin{\\left (x y \\right )}\\right)$$" ], "metadata": {}, "output_type": "pyout", "prompt_number": 54, "text": [ "-x\u22c5(x\u22c5y\u22c5cos(x\u22c5y) + 2\u22c5sin(x\u22c5y))" ] } ], "prompt_number": 54 }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Integration\n", "\n", "Integration is done in a similar fashion:" ] }, { "cell_type": "code", "collapsed": false, "input": [ "f" ], "language": "python", "metadata": {}, "outputs": [ { "latex": [ "$$\\sin{\\left (x y \\right )} + \\cos{\\left (y z \\right )}$$" ], "metadata": {}, "output_type": "pyout", "prompt_number": 55, "text": [ "sin(x\u22c5y) + cos(y\u22c5z)" ] } ], "prompt_number": 55 }, { "cell_type": "code", "collapsed": false, "input": [ "integrate(f, x)" ], "language": "python", "metadata": {}, "outputs": [ { "latex": [ "$$x \\cos{\\left (y z \\right )} + \\begin{cases} 0 & \\text{for}\\: y = 0 \\\\- \\frac{\\cos{\\left (x y \\right )}}{y} & \\text{otherwise} \\end{cases}$$" ], "metadata": {}, "output_type": "pyout", "prompt_number": 56, "text": [ " \u239b\u23a7 0 for y = 0\u239e\n", " \u239c\u23aa \u239f\n", "x\u22c5cos(y\u22c5z) + \u239c\u23a8-cos(x\u22c5y) \u239f\n", " \u239c\u23aa\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500 otherwise\u239f\n", " \u239d\u23a9 y \u23a0" ] } ], "prompt_number": 56 }, { "cell_type": "markdown", "metadata": {}, "source": [ "By providing limits for the integration variable we can evaluate definite integrals:" ] }, { "cell_type": "code", "collapsed": false, "input": [ "integrate(f, (x, -1, 1))" ], "language": "python", "metadata": {}, "outputs": [ { "latex": [ "$$2 \\cos{\\left (y z \\right )}$$" ], "metadata": {}, "output_type": "pyout", "prompt_number": 57, "text": [ "2\u22c5cos(y\u22c5z)" ] } ], "prompt_number": 57 }, { "cell_type": "markdown", "metadata": {}, "source": [ "and also improper integrals" ] }, { "cell_type": "code", "collapsed": false, "input": [ "integrate(exp(-x**2), (x, -oo, oo))" ], "language": "python", "metadata": {}, "outputs": [ { "latex": [ "$$\\sqrt{\\pi}$$" ], "metadata": {}, "output_type": "pyout", "prompt_number": 58, "text": [ " ___\n", "\u2572\u2571 \u03c0 " ] } ], "prompt_number": 58 }, { "cell_type": "markdown", "metadata": {}, "source": [ "Remember, `oo` is the SymPy notation for inifinity." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Sums and products\n", "\n", "We can evaluate sums and products using the functions: 'Sum'" ] }, { "cell_type": "code", "collapsed": false, "input": [ "n = Symbol(\"n\")" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 59 }, { "cell_type": "code", "collapsed": false, "input": [ "Sum(1/n**2, (n, 1, 10))" ], "language": "python", "metadata": {}, "outputs": [ { "latex": [ "$$\\sum_{n=1}^{10} n^{-2}$$" ], "metadata": {}, "output_type": "pyout", "prompt_number": 60, "text": [ " 10 \n", " ____ \n", " \u2572 \n", " \u2572 1 \n", " \u2572 \u2500\u2500\n", " \u2571 2\n", " \u2571 n \n", " \u2571 \n", " \u203e\u203e\u203e\u203e \n", "n = 1 " ] } ], "prompt_number": 60 }, { "cell_type": "code", "collapsed": false, "input": [ "Sum(1/n**2, (n,1, 10)).evalf()" ], "language": "python", "metadata": {}, "outputs": [ { "latex": [ "$$1.54976773116654$$" ], "metadata": {}, "output_type": "pyout", "prompt_number": 61, "text": [ "1.54976773116654" ] } ], "prompt_number": 61 }, { "cell_type": "code", "collapsed": false, "input": [ "Sum(1/n**2, (n, 1, oo)).evalf()" ], "language": "python", "metadata": {}, "outputs": [ { "latex": [ "$$1.64493406684823$$" ], "metadata": {}, "output_type": "pyout", "prompt_number": 62, "text": [ "1.64493406684823" ] } ], "prompt_number": 62 }, { "cell_type": "markdown", "metadata": {}, "source": [ "Products work much the same way:" ] }, { "cell_type": "code", "collapsed": false, "input": [ "Product(n, (n, 1, 10)) # 10!" ], "language": "python", "metadata": {}, "outputs": [ { "latex": [ "$$\\prod_{n=1}^{10} n$$" ], "metadata": {}, "output_type": "pyout", "prompt_number": 63, "text": [ " 10 \n", "\u252c\u2500\u2500\u2500\u252c \n", "\u2502 \u2502 n\n", "\u2502 \u2502 \n", "n = 1 " ] } ], "prompt_number": 63 }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Limits\n", "\n", "Limits can be evaluated using the `limit` function. For example, " ] }, { "cell_type": "code", "collapsed": false, "input": [ "limit(sin(x)/x, x, 0)" ], "language": "python", "metadata": {}, "outputs": [ { "latex": [ "$$1$$" ], "metadata": {}, "output_type": "pyout", "prompt_number": 64, "text": [ "1" ] } ], "prompt_number": 64 }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can use 'limit' to check the result of derivation using the `diff` function:" ] }, { "cell_type": "code", "collapsed": false, "input": [ "f" ], "language": "python", "metadata": {}, "outputs": [ { "latex": [ "$$\\sin{\\left (x y \\right )} + \\cos{\\left (y z \\right )}$$" ], "metadata": {}, "output_type": "pyout", "prompt_number": 65, "text": [ "sin(x\u22c5y) + cos(y\u22c5z)" ] } ], "prompt_number": 65 }, { "cell_type": "code", "collapsed": false, "input": [ "diff(f, x)" ], "language": "python", "metadata": {}, "outputs": [ { "latex": [ "$$y \\cos{\\left (x y \\right )}$$" ], "metadata": {}, "output_type": "pyout", "prompt_number": 66, "text": [ "y\u22c5cos(x\u22c5y)" ] } ], "prompt_number": 66 }, { "cell_type": "markdown", "metadata": {}, "source": [ "$\\displaystyle \\frac{\\mathrm{d}f(x,y)}{\\mathrm{d}x} = \\frac{f(x+h,y)-f(x,y)}{h}$" ] }, { "cell_type": "code", "collapsed": false, "input": [ "h = Symbol(\"h\")" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 67 }, { "cell_type": "code", "collapsed": false, "input": [ "limit((f.subs(x, x+h) - f)/h, h, 0)" ], "language": "python", "metadata": {}, "outputs": [ { "latex": [ "$$y \\cos{\\left (x y \\right )}$$" ], "metadata": {}, "output_type": "pyout", "prompt_number": 68, "text": [ "y\u22c5cos(x\u22c5y)" ] } ], "prompt_number": 68 }, { "cell_type": "markdown", "metadata": {}, "source": [ "OK!" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can change the direction from which we approach the limiting point using the `dir` keywork argument:" ] }, { "cell_type": "code", "collapsed": false, "input": [ "limit(1/x, x, 0, dir=\"+\")" ], "language": "python", "metadata": {}, "outputs": [ { "latex": [ "$$\\infty$$" ], "metadata": {}, "output_type": "pyout", "prompt_number": 69, "text": [ "\u221e" ] } ], "prompt_number": 69 }, { "cell_type": "code", "collapsed": false, "input": [ "limit(1/x, x, 0, dir=\"-\")" ], "language": "python", "metadata": {}, "outputs": [ { "latex": [ "$$-\\infty$$" ], "metadata": {}, "output_type": "pyout", "prompt_number": 70, "text": [ "-\u221e" ] } ], "prompt_number": 70 }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Series\n", "\n", "Series expansion is also one of the most useful features of a CAS. In SymPy we can perform a series expansion of an expression using the `series` function:" ] }, { "cell_type": "code", "collapsed": false, "input": [ "series(exp(x), x)" ], "language": "python", "metadata": {}, "outputs": [ { "latex": [ "$$1 + x + \\frac{1}{2} x^{2} + \\frac{1}{6} x^{3} + \\frac{1}{24} x^{4} + \\frac{1}{120} x^{5} + \\mathcal{O}\\left(x^{6}\\right)$$" ], "metadata": {}, "output_type": "pyout", "prompt_number": 71, "text": [ " 2 3 4 5 \n", " x x x x \u239b 6\u239e\n", "1 + x + \u2500\u2500 + \u2500\u2500 + \u2500\u2500 + \u2500\u2500\u2500 + O\u239dx \u23a0\n", " 2 6 24 120 " ] } ], "prompt_number": 71 }, { "cell_type": "markdown", "metadata": {}, "source": [ "By default it expands the expression around $x=0$, but we can expand around any value of $x$ by explicitly include a value in the function call:" ] }, { "cell_type": "code", "collapsed": false, "input": [ "series(exp(x), x, 1)" ], "language": "python", "metadata": {}, "outputs": [ { "latex": [ "$$e + e x + \\frac{1}{2} e x^{2} + \\frac{1}{6} e x^{3} + \\frac{1}{24} e x^{4} + \\frac{1}{120} e x^{5} + \\mathcal{O}\\left(x^{6}\\right)$$" ], "metadata": {}, "output_type": "pyout", "prompt_number": 72, "text": [ " 2 3 4 5 \n", " \u212f\u22c5x \u212f\u22c5x \u212f\u22c5x \u212f\u22c5x \u239b 6\u239e\n", "\u212f + \u212f\u22c5x + \u2500\u2500\u2500\u2500 + \u2500\u2500\u2500\u2500 + \u2500\u2500\u2500\u2500 + \u2500\u2500\u2500\u2500 + O\u239dx \u23a0\n", " 2 6 24 120 " ] } ], "prompt_number": 72 }, { "cell_type": "markdown", "metadata": {}, "source": [ "And we can explicitly define to which order the series expansion should be carried out:" ] }, { "cell_type": "code", "collapsed": false, "input": [ "series(exp(x), x, 1, 10)" ], "language": "python", "metadata": {}, "outputs": [ { "latex": [ "$$e + e x + \\frac{1}{2} e x^{2} + \\frac{1}{6} e x^{3} + \\frac{1}{24} e x^{4} + \\frac{1}{120} e x^{5} + \\frac{1}{720} e x^{6} + \\frac{1}{5040} e x^{7} + \\frac{1}{40320} e x^{8} + \\frac{1}{362880} e x^{9} + \\mathcal{O}\\left(x^{10}\\right)$$" ], "metadata": {}, "output_type": "pyout", "prompt_number": 73, "text": [ " 2 3 4 5 6 7 8 9 \n", " \u212f\u22c5x \u212f\u22c5x \u212f\u22c5x \u212f\u22c5x \u212f\u22c5x \u212f\u22c5x \u212f\u22c5x \u212f\u22c5x \u239b 10\u239e\n", "\u212f + \u212f\u22c5x + \u2500\u2500\u2500\u2500 + \u2500\u2500\u2500\u2500 + \u2500\u2500\u2500\u2500 + \u2500\u2500\u2500\u2500 + \u2500\u2500\u2500\u2500 + \u2500\u2500\u2500\u2500 + \u2500\u2500\u2500\u2500\u2500 + \u2500\u2500\u2500\u2500\u2500\u2500 + O\u239dx \u23a0\n", " 2 6 24 120 720 5040 40320 362880 " ] } ], "prompt_number": 73 }, { "cell_type": "markdown", "metadata": {}, "source": [ "The series expansion includes the order of the approximation, which is very useful for keeping track of the order of validity when we do calculations with series expansions of different order:" ] }, { "cell_type": "code", "collapsed": false, "input": [ "s1 = cos(x).series(x, 0, 5)\n", "s1" ], "language": "python", "metadata": {}, "outputs": [ { "latex": [ "$$1 - \\frac{1}{2} x^{2} + \\frac{1}{24} x^{4} + \\mathcal{O}\\left(x^{5}\\right)$$" ], "metadata": {}, "output_type": "pyout", "prompt_number": 74, "text": [ " 2 4 \n", " x x \u239b 5\u239e\n", "1 - \u2500\u2500 + \u2500\u2500 + O\u239dx \u23a0\n", " 2 24 " ] } ], "prompt_number": 74 }, { "cell_type": "code", "collapsed": false, "input": [ "s2 = sin(x).series(x, 0, 2)\n", "s2" ], "language": "python", "metadata": {}, "outputs": [ { "latex": [ "$$x + \\mathcal{O}\\left(x^{2}\\right)$$" ], "metadata": {}, "output_type": "pyout", "prompt_number": 75, "text": [ " \u239b 2\u239e\n", "x + O\u239dx \u23a0" ] } ], "prompt_number": 75 }, { "cell_type": "code", "collapsed": false, "input": [ "expand(s1 * s2)" ], "language": "python", "metadata": {}, "outputs": [ { "latex": [ "$$x + \\mathcal{O}\\left(x^{2}\\right)$$" ], "metadata": {}, "output_type": "pyout", "prompt_number": 76, "text": [ " \u239b 2\u239e\n", "x + O\u239dx \u23a0" ] } ], "prompt_number": 76 }, { "cell_type": "markdown", "metadata": {}, "source": [ "If we want to get rid of the order information we can use the `removeO` method:" ] }, { "cell_type": "code", "collapsed": false, "input": [ "expand(s1.removeO() * s2.removeO())" ], "language": "python", "metadata": {}, "outputs": [ { "latex": [ "$$\\frac{1}{24} x^{5} - \\frac{1}{2} x^{3} + x$$" ], "metadata": {}, "output_type": "pyout", "prompt_number": 77, "text": [ " 5 3 \n", "x x \n", "\u2500\u2500 - \u2500\u2500 + x\n", "24 2 " ] } ], "prompt_number": 77 }, { "cell_type": "markdown", "metadata": {}, "source": [ "But note that this is not the correct expansion of $\\cos(x)\\sin(x)$ to $5$th order:" ] }, { "cell_type": "code", "collapsed": false, "input": [ "(cos(x)*sin(x)).series(x, 0, 6)" ], "language": "python", "metadata": {}, "outputs": [ { "latex": [ "$$x - \\frac{2}{3} x^{3} + \\frac{2}{15} x^{5} + \\mathcal{O}\\left(x^{6}\\right)$$" ], "metadata": {}, "output_type": "pyout", "prompt_number": 78, "text": [ " 3 5 \n", " 2\u22c5x 2\u22c5x \u239b 6\u239e\n", "x - \u2500\u2500\u2500\u2500 + \u2500\u2500\u2500\u2500 + O\u239dx \u23a0\n", " 3 15 " ] } ], "prompt_number": 78 }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Linear algebra" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Matrices\n", "\n", "Matrices are defined using the `Matrix` class:" ] }, { "cell_type": "code", "collapsed": false, "input": [ "m11, m12, m21, m22 = symbols(\"m11, m12, m21, m22\")\n", "b1, b2 = symbols(\"b1, b2\")" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 79 }, { "cell_type": "code", "collapsed": false, "input": [ "A = Matrix([[m11, m12],[m21, m22]])\n", "A" ], "language": "python", "metadata": {}, "outputs": [ { "latex": [ "$$\\left[\\begin{smallmatrix}m_{11} & m_{12}\\\\m_{21} & m_{22}\\end{smallmatrix}\\right]$$" ], "metadata": {}, "output_type": "pyout", "prompt_number": 80, "text": [ "\u23a1m\u2081\u2081 m\u2081\u2082\u23a4\n", "\u23a2 \u23a5\n", "\u23a3m\u2082\u2081 m\u2082\u2082\u23a6" ] } ], "prompt_number": 80 }, { "cell_type": "code", "collapsed": false, "input": [ "b = Matrix([[b1], [b2]])\n", "b" ], "language": "python", "metadata": {}, "outputs": [ { "latex": [ "$$\\left[\\begin{smallmatrix}b_{1}\\\\b_{2}\\end{smallmatrix}\\right]$$" ], "metadata": {}, "output_type": "pyout", "prompt_number": 81, "text": [ "\u23a1b\u2081\u23a4\n", "\u23a2 \u23a5\n", "\u23a3b\u2082\u23a6" ] } ], "prompt_number": 81 }, { "cell_type": "markdown", "metadata": {}, "source": [ "With `Matrix` class instances we can do the usual matrix algebra operations:" ] }, { "cell_type": "code", "collapsed": false, "input": [ "A**2" ], "language": "python", "metadata": {}, "outputs": [ { "latex": [ "$$\\left[\\begin{smallmatrix}m_{11}^{2} + m_{12} m_{21} & m_{11} m_{12} + m_{12} m_{22}\\\\m_{11} m_{21} + m_{21} m_{22} & m_{12} m_{21} + m_{22}^{2}\\end{smallmatrix}\\right]$$" ], "metadata": {}, "output_type": "pyout", "prompt_number": 82, "text": [ "\u23a1 2 \u23a4\n", "\u23a2 m\u2081\u2081 + m\u2081\u2082\u22c5m\u2082\u2081 m\u2081\u2081\u22c5m\u2081\u2082 + m\u2081\u2082\u22c5m\u2082\u2082\u23a5\n", "\u23a2 \u23a5\n", "\u23a2 2 \u23a5\n", "\u23a3m\u2081\u2081\u22c5m\u2082\u2081 + m\u2082\u2081\u22c5m\u2082\u2082 m\u2081\u2082\u22c5m\u2082\u2081 + m\u2082\u2082 \u23a6" ] } ], "prompt_number": 82 }, { "cell_type": "code", "collapsed": false, "input": [ "A * b" ], "language": "python", "metadata": {}, "outputs": [ { "latex": [ "$$\\left[\\begin{smallmatrix}b_{1} m_{11} + b_{2} m_{12}\\\\b_{1} m_{21} + b_{2} m_{22}\\end{smallmatrix}\\right]$$" ], "metadata": {}, "output_type": "pyout", "prompt_number": 83, "text": [ "\u23a1b\u2081\u22c5m\u2081\u2081 + b\u2082\u22c5m\u2081\u2082\u23a4\n", "\u23a2 \u23a5\n", "\u23a3b\u2081\u22c5m\u2082\u2081 + b\u2082\u22c5m\u2082\u2082\u23a6" ] } ], "prompt_number": 83 }, { "cell_type": "markdown", "metadata": {}, "source": [ "And calculate determinants and inverses, and the like:" ] }, { "cell_type": "code", "collapsed": false, "input": [ "A.det()" ], "language": "python", "metadata": {}, "outputs": [ { "latex": [ "$$m_{11} m_{22} - m_{12} m_{21}$$" ], "metadata": {}, "output_type": "pyout", "prompt_number": 84, "text": [ "m\u2081\u2081\u22c5m\u2082\u2082 - m\u2081\u2082\u22c5m\u2082\u2081" ] } ], "prompt_number": 84 }, { "cell_type": "code", "collapsed": false, "input": [ "A.inv()" ], "language": "python", "metadata": {}, "outputs": [ { "latex": [ "$$\\left[\\begin{smallmatrix}\\frac{1}{m_{11}} + \\frac{m_{12} m_{21}}{m_{11}^{2} \\left(m_{22} - \\frac{m_{12} m_{21}}{m_{11}}\\right)} & - \\frac{m_{12}}{m_{11} \\left(m_{22} - \\frac{m_{12} m_{21}}{m_{11}}\\right)}\\\\- \\frac{m_{21}}{m_{11} \\left(m_{22} - \\frac{m_{12} m_{21}}{m_{11}}\\right)} & \\frac{1}{m_{22} - \\frac{m_{12} m_{21}}{m_{11}}}\\end{smallmatrix}\\right]$$" ], "metadata": {}, "output_type": "pyout", "prompt_number": 85, "text": [ "\u23a1 1 m\u2081\u2082\u22c5m\u2082\u2081 -m\u2081\u2082 \u23a4\n", "\u23a2\u2500\u2500\u2500 + \u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500 \u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u23a5\n", "\u23a2m\u2081\u2081 2 \u239b m\u2081\u2082\u22c5m\u2082\u2081\u239e \u239b m\u2081\u2082\u22c5m\u2082\u2081\u239e\u23a5\n", "\u23a2 m\u2081\u2081 \u22c5\u239cm\u2082\u2082 - \u2500\u2500\u2500\u2500\u2500\u2500\u2500\u239f m\u2081\u2081\u22c5\u239cm\u2082\u2082 - \u2500\u2500\u2500\u2500\u2500\u2500\u2500\u239f\u23a5\n", "\u23a2 \u239d m\u2081\u2081 \u23a0 \u239d m\u2081\u2081 \u23a0\u23a5\n", "\u23a2 \u23a5\n", "\u23a2 -m\u2082\u2081 1 \u23a5\n", "\u23a2 \u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500 \u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500 \u23a5\n", "\u23a2 \u239b m\u2081\u2082\u22c5m\u2082\u2081\u239e m\u2081\u2082\u22c5m\u2082\u2081 \u23a5\n", "\u23a2 m\u2081\u2081\u22c5\u239cm\u2082\u2082 - \u2500\u2500\u2500\u2500\u2500\u2500\u2500\u239f m\u2082\u2082 - \u2500\u2500\u2500\u2500\u2500\u2500\u2500 \u23a5\n", "\u23a3 \u239d m\u2081\u2081 \u23a0 m\u2081\u2081 \u23a6" ] } ], "prompt_number": 85 }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Solving equations\n", "\n", "For solving equations and systems of equations we can use the `solve` function:" ] }, { "cell_type": "code", "collapsed": false, "input": [ "solve(x**2 - 1, x)" ], "language": "python", "metadata": {}, "outputs": [ { "latex": [ "$$\\begin{bmatrix}-1, & 1\\end{bmatrix}$$" ], "metadata": {}, "output_type": "pyout", "prompt_number": 86, "text": [ "[-1, 1]" ] } ], "prompt_number": 86 }, { "cell_type": "code", "collapsed": false, "input": [ "solve(x**4 - x**2 - 1, x)" ], "language": "python", "metadata": {}, "outputs": [ { "latex": [ "$$\\begin{bmatrix}- i \\sqrt{- \\frac{1}{2} + \\frac{1}{2} \\sqrt{5}}, & i \\sqrt{- \\frac{1}{2} + \\frac{1}{2} \\sqrt{5}}, & - \\sqrt{\\frac{1}{2} + \\frac{1}{2} \\sqrt{5}}, & \\sqrt{\\frac{1}{2} + \\frac{1}{2} \\sqrt{5}}\\end{bmatrix}$$" ], "metadata": {}, "output_type": "pyout", "prompt_number": 87, "text": [ "\u23a1 _____________ _____________ ___________ ________\n", "\u23a2 \u2571 ___ \u2571 ___ \u2571 ___ \u2571 _\n", "\u23a2 \u2571 1 \u2572\u2571 5 \u2571 1 \u2572\u2571 5 \u2571 1 \u2572\u2571 5 \u2571 1 \u2572\u2571 \n", "\u23a2-\u2148\u22c5 \u2571 - \u2500 + \u2500\u2500\u2500\u2500\u2500 , \u2148\u22c5 \u2571 - \u2500 + \u2500\u2500\u2500\u2500\u2500 , - \u2571 \u2500 + \u2500\u2500\u2500\u2500\u2500 , \u2571 \u2500 + \u2500\u2500\u2500\n", "\u23a3 \u2572\u2571 2 2 \u2572\u2571 2 2 \u2572\u2571 2 2 \u2572\u2571 2 2\n", "\n", "___\u23a4\n", "__ \u23a5\n", "5 \u23a5\n", "\u2500\u2500 \u23a5\n", " \u23a6" ] } ], "prompt_number": 87 }, { "cell_type": "markdown", "metadata": {}, "source": [ "System of equations:" ] }, { "cell_type": "code", "collapsed": false, "input": [ "solve([x + y - 1, x - y - 1], [x,y])" ], "language": "python", "metadata": {}, "outputs": [ { "latex": [ "$$\\begin{Bmatrix}x : 1, & y : 0\\end{Bmatrix}$$" ], "metadata": {}, "output_type": "pyout", "prompt_number": 88, "text": [ "{x: 1, y: 0}" ] } ], "prompt_number": 88 }, { "cell_type": "markdown", "metadata": {}, "source": [ "In terms of other symbolic expressions:" ] }, { "cell_type": "code", "collapsed": false, "input": [ "solve([x + y - a, x - y - c], [x,y])" ], "language": "python", "metadata": {}, "outputs": [ { "latex": [ "$$\\begin{Bmatrix}x : \\frac{1}{2} a + \\frac{1}{2} c, & y : \\frac{1}{2} a - \\frac{1}{2} c\\end{Bmatrix}$$" ], "metadata": {}, "output_type": "pyout", "prompt_number": 89, "text": [ "\u23a7 a c a c\u23ab\n", "\u23a8x: \u2500 + \u2500, y: \u2500 - \u2500\u23ac\n", "\u23a9 2 2 2 2\u23ad" ] } ], "prompt_number": 89 }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Quantum mechanics: noncommuting variables\n", "\n", "How about non-commuting symbols? In quantum mechanics we need to work with noncommuting operators, and SymPy has a nice support for noncommuting symbols and even a subpackage for quantum mechanics related calculations!" ] }, { "cell_type": "code", "collapsed": false, "input": [ "from sympy.physics.quantum import *" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 90 }, { "cell_type": "markdown", "metadata": {}, "source": [ "### States\n", "\n", "We can define symbol states, kets and bras:" ] }, { "cell_type": "code", "collapsed": false, "input": [ "Ket('psi')" ], "language": "python", "metadata": {}, "outputs": [ { "latex": [ "$${\\left|\\psi\\right\\rangle }$$" ], "metadata": {}, "output_type": "pyout", "prompt_number": 91, "text": [ "\u2758\u03c8\u27e9" ] } ], "prompt_number": 91 }, { "cell_type": "code", "collapsed": false, "input": [ "Bra('psi')" ], "language": "python", "metadata": {}, "outputs": [ { "latex": [ "$${\\left\\langle \\psi\\right|}$$" ], "metadata": {}, "output_type": "pyout", "prompt_number": 92, "text": [ "\u27e8\u03c8\u2758" ] } ], "prompt_number": 92 }, { "cell_type": "code", "collapsed": false, "input": [ "u = Ket('0')\n", "d = Ket('1')\n", "\n", "a, b = symbols('alpha beta', complex=True)" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 93 }, { "cell_type": "code", "collapsed": false, "input": [ "phi = a * u + sqrt(1-abs(a)**2) * d; phi" ], "language": "python", "metadata": {}, "outputs": [ { "latex": [ "$$\\alpha {\\left|0\\right\\rangle } + \\sqrt{- \\left\\lvert{\\alpha}\\right\\rvert^{2} + 1} {\\left|1\\right\\rangle }$$" ], "metadata": {}, "output_type": "pyout", "prompt_number": 94, "text": [ " ____________ \n", " \u2571 2 \n", "\u03b1\u22c5\u27580\u27e9 + \u2572\u2571 - \u2502\u03b1\u2502 + 1 \u22c5\u27581\u27e9" ] } ], "prompt_number": 94 }, { "cell_type": "code", "collapsed": false, "input": [ "Dagger(phi)" ], "language": "python", "metadata": {}, "outputs": [ { "latex": [ "$$\\overline{\\alpha} {\\left\\langle 0\\right|} + \\overline{\\sqrt{- \\left\\lvert{\\alpha}\\right\\rvert^{2} + 1}} {\\left\\langle 1\\right|}$$" ], "metadata": {}, "output_type": "pyout", "prompt_number": 95, "text": [ " _______________ \n", " ____________ \n", "_ \u2571 2 \n", "\u03b1\u22c5\u27e80\u2758 + \u2572\u2571 - \u2502\u03b1\u2502 + 1 \u22c5\u27e81\u2758" ] } ], "prompt_number": 95 }, { "cell_type": "code", "collapsed": false, "input": [ "Dagger(phi) * d" ], "language": "python", "metadata": {}, "outputs": [ { "latex": [ "$$\\left(\\overline{\\alpha} {\\left\\langle 0\\right|} + \\overline{\\sqrt{- \\left\\lvert{\\alpha}\\right\\rvert^{2} + 1}} {\\left\\langle 1\\right|}\\right) {\\left|1\\right\\rangle }$$" ], "metadata": {}, "output_type": "pyout", "prompt_number": 96, "text": [ "\u239b _______________ \u239e \n", "\u239c ____________ \u239f \n", "\u239c_ \u2571 2 \u239f \n", "\u239d\u03b1\u22c5\u27e80\u2758 + \u2572\u2571 - \u2502\u03b1\u2502 + 1 \u22c5\u27e81\u2758\u23a0\u22c5\u27581\u27e9" ] } ], "prompt_number": 96 }, { "cell_type": "markdown", "metadata": {}, "source": [ "Use `qapply` to distribute a mutiplication:" ] }, { "cell_type": "code", "collapsed": false, "input": [ "qapply(Dagger(phi) * d)" ], "language": "python", "metadata": {}, "outputs": [ { "latex": [ "$$\\overline{\\alpha} \\left\\langle 0 \\right. {\\left|1\\right\\rangle } + \\overline{\\sqrt{- \\left\\lvert{\\alpha}\\right\\rvert^{2} + 1}} \\left\\langle 1 \\right. {\\left|1\\right\\rangle }$$" ], "metadata": {}, "output_type": "pyout", "prompt_number": 97, "text": [ " _______________ \n", " ____________ \n", "_ \u2571 2 \n", "\u03b1\u22c5\u27e80\u27581\u27e9 + \u2572\u2571 - \u2502\u03b1\u2502 + 1 \u22c5\u27e81\u27581\u27e9" ] } ], "prompt_number": 97 }, { "cell_type": "code", "collapsed": false, "input": [ "qapply(Dagger(phi) * u)" ], "language": "python", "metadata": {}, "outputs": [ { "latex": [ "$$\\overline{\\alpha} \\left\\langle 0 \\right. {\\left|0\\right\\rangle } + \\overline{\\sqrt{- \\left\\lvert{\\alpha}\\right\\rvert^{2} + 1}} \\left\\langle 1 \\right. {\\left|0\\right\\rangle }$$" ], "metadata": {}, "output_type": "pyout", "prompt_number": 98, "text": [ " _______________ \n", " ____________ \n", "_ \u2571 2 \n", "\u03b1\u22c5\u27e80\u27580\u27e9 + \u2572\u2571 - \u2502\u03b1\u2502 + 1 \u22c5\u27e81\u27580\u27e9" ] } ], "prompt_number": 98 }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Operators" ] }, { "cell_type": "code", "collapsed": false, "input": [ "A = Operator('A')\n", "B = Operator('B')" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 99 }, { "cell_type": "markdown", "metadata": {}, "source": [ "Check if they are commuting!" ] }, { "cell_type": "code", "collapsed": false, "input": [ "A * B == B * A" ], "language": "python", "metadata": {}, "outputs": [ { "latex": [ "$$False$$" ], "metadata": {}, "output_type": "pyout", "prompt_number": 100, "text": [ "False" ] } ], "prompt_number": 100 }, { "cell_type": "code", "collapsed": false, "input": [ "expand((A+B)**3)" ], "language": "python", "metadata": {}, "outputs": [ { "latex": [ "$$A B A + A \\left(B\\right)^{2} + \\left(A\\right)^{2} B + \\left(A\\right)^{3} + B A B + B \\left(A\\right)^{2} + \\left(B\\right)^{2} A + \\left(B\\right)^{3}$$" ], "metadata": {}, "output_type": "pyout", "prompt_number": 101, "text": [ " 2 2 3 2 2 3\n", "A\u22c5B\u22c5A + A\u22c5B + A \u22c5B + A + B\u22c5A\u22c5B + B\u22c5A + B \u22c5A + B " ] } ], "prompt_number": 101 }, { "cell_type": "code", "collapsed": false, "input": [ "c = Commutator(A,B)\n", "c" ], "language": "python", "metadata": {}, "outputs": [ { "latex": [ "$$\\left[A,B\\right]$$" ], "metadata": {}, "output_type": "pyout", "prompt_number": 102, "text": [ "[A,B]" ] } ], "prompt_number": 102 }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can use the `doit` method to evaluate the commutator:" ] }, { "cell_type": "code", "collapsed": false, "input": [ "c.doit()" ], "language": "python", "metadata": {}, "outputs": [ { "latex": [ "$$A B - B A$$" ], "metadata": {}, "output_type": "pyout", "prompt_number": 103, "text": [ "A\u22c5B - B\u22c5A" ] } ], "prompt_number": 103 }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can mix quantum operators with C-numbers:" ] }, { "cell_type": "code", "collapsed": false, "input": [ "c = Commutator(a * A, b * B)\n", "c" ], "language": "python", "metadata": {}, "outputs": [ { "latex": [ "$$\\alpha \\beta \\left[A,B\\right]$$" ], "metadata": {}, "output_type": "pyout", "prompt_number": 104, "text": [ "\u03b1\u22c5\u03b2\u22c5[A,B]" ] } ], "prompt_number": 104 }, { "cell_type": "markdown", "metadata": {}, "source": [ "To expand the commutator, use the `expand` method with the `commutator=True` keyword argument:" ] }, { "cell_type": "code", "collapsed": false, "input": [ "c = Commutator(A+B, A*B)\n", "c.expand(commutator=True)" ], "language": "python", "metadata": {}, "outputs": [ { "latex": [ "$$- \\left[A,B\\right] B + A \\left[A,B\\right]$$" ], "metadata": {}, "output_type": "pyout", "prompt_number": 105, "text": [ "-[A,B]\u22c5B + A\u22c5[A,B]" ] } ], "prompt_number": 105 }, { "cell_type": "code", "collapsed": false, "input": [ "Dagger(Commutator(A, B))" ], "language": "python", "metadata": {}, "outputs": [ { "latex": [ "$$- \\left[A^{\\dagger},B^{\\dagger}\\right]$$" ], "metadata": {}, "output_type": "pyout", "prompt_number": 106, "text": [ " \u23a1 \u2020 \u2020\u23a4\n", "-\u23a3A ,B \u23a6" ] } ], "prompt_number": 106 }, { "cell_type": "code", "collapsed": false, "input": [ "ac = AntiCommutator(A,B)" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 107 }, { "cell_type": "code", "collapsed": false, "input": [ "ac.doit()" ], "language": "python", "metadata": {}, "outputs": [ { "latex": [ "$$A B + B A$$" ], "metadata": {}, "output_type": "pyout", "prompt_number": 108, "text": [ "A\u22c5B + B\u22c5A" ] } ], "prompt_number": 108 }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Example: Quadrature commutator" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let's look at the commutator of the electromagnetic field quadatures $x$ and $p$. We can write the quadrature operators in terms of the creation and annihilation operators as:\n", "\n", "$\\displaystyle x = (a + a^\\dagger)/\\sqrt{2}$\n", "\n", "$\\displaystyle p = -i(a - a^\\dagger)/\\sqrt{2}$\n" ] }, { "cell_type": "code", "collapsed": false, "input": [ "X = (A + Dagger(A))/sqrt(2)\n", "X" ], "language": "python", "metadata": {}, "outputs": [ { "latex": [ "$$\\frac{1}{2} \\sqrt{2} \\left(A^{\\dagger} + A\\right)$$" ], "metadata": {}, "output_type": "pyout", "prompt_number": 109, "text": [ " ___ \u239b \u2020 \u239e\n", "\u2572\u2571 2 \u22c5\u239dA + A\u23a0\n", "\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\n", " 2 " ] } ], "prompt_number": 109 }, { "cell_type": "code", "collapsed": false, "input": [ "P = -I * (A - Dagger(A))/sqrt(2)\n", "P" ], "language": "python", "metadata": {}, "outputs": [ { "latex": [ "$$- \\frac{1}{2} \\sqrt{2} i \\left(- A^{\\dagger} + A\\right)$$" ], "metadata": {}, "output_type": "pyout", "prompt_number": 110, "text": [ " ___ \u239b \u2020 \u239e\n", "-\u2572\u2571 2 \u22c5\u2148\u22c5\u239d- A + A\u23a0\n", "\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\n", " 2 " ] } ], "prompt_number": 110 }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let's expand the commutator $[x,p]$" ] }, { "cell_type": "code", "collapsed": false, "input": [ "Commutator(X, P).expand(commutator=True).expand(commutator=True)" ], "language": "python", "metadata": {}, "outputs": [ { "latex": [ "$$- i \\left[A^{\\dagger},A\\right]$$" ], "metadata": {}, "output_type": "pyout", "prompt_number": 111, "text": [ " \u23a1 \u2020 \u23a4\n", "-\u2148\u22c5\u23a3A ,A\u23a6" ] } ], "prompt_number": 111 }, { "cell_type": "markdown", "metadata": {}, "source": [ "Here we see directly that the well known commutation relation for the quadratures\n", "\n", "$[x,p]=i$\n", "\n", "is a directly related to\n", "\n", "$[A, A^\\dagger]=1$ \n", "\n", "(which SymPy does not know about, and does not simplify)." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "For more details on the quantum module in SymPy, see:\n", "\n", "* http://docs.sympy.org/0.7.2/modules/physics/quantum/index.html\n", "* http://nbviewer.ipython.org/urls/raw.github.com/ipython/ipython/master/docs/examples/notebooks/sympy_quantum_computing.ipynb" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Further reading\n", "\n", "* http://sympy.org/en/index.html - The SymPy projects web page.\n", "* https://github.com/sympy/sympy - The source code of SymPy.\n", "* http://live.sympy.org - Online version of SymPy for testing and demonstrations.\n" ] } ], "metadata": {} } ] }