{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Symbolic computation with SymPy\n", "In the previous chapters we have studied methods for solving problems numerically, i.e. to find approximate solutions in particular when the analytical solution is impossible or tedious. \n", "\n", "The [SymPy package](https://docs.sympy.org/) (a Python library for symbolic mathematics) offers an alternative approach and allows to do symbolic computation, i.e. the algebraic manipulation of expressions, in Python in a very natural way. " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Basics" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let us start with a very simple example that demonstrates the difference between numerical and symbolic computation. To compute the square root of a number, we can use the `sqrt` function from the `math` module: " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import math\n", "math.sqrt(8)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The result is an irrational number, which we know is not only $2.828...$ but can also be written as $\\sqrt8 = \\sqrt{4\\cdot2} = 2\\cdot\\sqrt2$. However, from the result of the `math.sqrt` function, $2.828...$, this last relation is not obvious. Wouldn't it be nice if Python could tell us this?" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import sympy\n", "sympy.init_printing() # enable LaTeX output (or best available)\n", "sympy.sqrt(8)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "This is not limited to numerical expressions but can also be used for algebraic expressions. \n", "For these, we need to define the symbols we want to use:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "from sympy import symbols\n", "x, y = symbols('x y')\n", "x" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Not only is this symbol now rendered as LaTeX; it's no longer a simple Python variable." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "type(x)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "(Question: what would this output \"normally\", i.e. after `x = 5`?)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can now use it in expressions without assigning a specific value:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "f = x + 2*y\n", "f" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "(Question: what would this output if we had not defined `x` before?)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "f - x" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "fsquare = f**2\n", "fsquare" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "To expand the brackets in this symbolic expression, we can use `expand`:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "from sympy import expand\n", "g = expand(fsquare)\n", "g" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The inverse of `expand` is `simplify` or, for this specific case, `factor`:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "from sympy import simplify, factor\n", "factor(g)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Note that `fsquare`$ = (x+2y)^2$ and `g`$ = x^2+4xy+4y^2$ are [not](https://docs.sympy.org/dev/tutorials/intro-tutorial/gotchas.html#equals-signs) considered to be identical:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "fsquare == g" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "However, we can compute their difference and see that this is zero:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "fsquare - g" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now, is this zero? 🤔\n", "\n", "Let us simplify it to see:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "simplify(fsquare - g)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "There is an alternative way:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "fsquare.equals(g)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Note though that this alternative does not test equality symbolically but numerically at random points." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "To evaluate an expression with a concrete value, we can use `subs` (for substitute):" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "f.subs(x, 3)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "or for several values" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "f.subs([(x, 3), (y, 4)])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "This evaluates to a simple number as you would expect. What about the $\\sqrt8$ from above? For this, we can also compute the numeric value using `evalf`:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "from sympy import sqrt, evalf\n", "w8 = sqrt(8)\n", "w8" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "w8.evalf() # optional: specify accuracy as argument" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "To use numbers as algebraic expressions, e.g. in fractions, we can use `S()`:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "from sympy import S\n", "print(f\"1/2 = {1/2}\")\n", "print(f\"{S(1)/2 = }\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "(Question: would `S(1/2)` work as well?)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "SymPy also provides standard math functions of which we're going to use a couple in the following:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "from sympy import exp, sin, cos, tan, sinh" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Calculus and Linear Algebra in SymPy " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "So far, we have only scratched the surface with the basic algebraic transformations above, and there is no time here to discuss all functionality (and pitfalls) in depth, so we will limit ourselves to some of the most important further applications here." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Differentiation\n", "SymPy can be used to find derivatives:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "from sympy import diff\n", "diff(x**2) # implies derivative w.r.t. x" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Multiple times:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "diff(x**2, x, x) # alternative: diff(x**2, x, 2)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Multiple variables:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "f = x**2+x*sin(y)\n", "f" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Derivative w.r.t. $y$:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "diff(f, y)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Integration\n", "SymPy can be used to find integrals:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "from sympy import integrate\n", "integrate(2*x)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Indefinite integral of the standard normal distribution, $\\int_{-\\infty}^{\\infty} \\frac{\\exp(-x^2/2)}{\\sqrt{2\\pi}}$:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "gauss = exp(-x**2 / 2) / sqrt(2*sympy.pi)\n", "integrate(gauss)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Or its integral over the real numbers with the symbol `oo` for $\\infty$:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "from sympy import oo\n", "integrate(gauss, (x, -oo, oo))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Taking limits\n", "Compute limits, e.g. $\\lim_{x\\to0} \\frac{\\sin(x)}x$:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "from sympy import limit\n", "\n", "limit(sin(x) / x, x, 0)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Rewriting\n", "Simplifying functions by rewriting:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "f = exp(x) * sinh(x)\n", "f" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "f_exp = f.rewrite(exp)\n", "f_exp" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "f_exp.simplify()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Direct way does not work:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "f.simplify()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Also useful to rewrite trigonometrical expressions:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "f = sin(2*x)*tan(x**2)\n", "f" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "f.rewrite(sin)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Series expansions" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "from sympy import series\n", "sin(x).series(x, 0, 10)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Solving equations\n", "Symbolic relations in SymPy are not written with `=` (assignment) or `==` (comparison) but with the special keyword `Eq`:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "from sympy import Eq\n", "Eq(x, y)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can use this to solve the equation $x^2 = -1$ (over the complex or real numbers) in SymPy:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "from sympy import solveset\n", "solveset(Eq(x**2, -1), x)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "solveset(Eq(x**2, -1), x, domain = S.Reals)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "There is a equivalent short-hand notation to solve \"expression = 0\":" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "solveset(x**2 + 1, x)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Linear algebra\n", "Compute eigenvalues of a matrix:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "from sympy import Matrix\n", "M = Matrix([[1, 2], [2, 1]])\n", "M" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "l1, l2 = M.eigenvals()\n", "l1, l2" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "M.eigenvects() # eigenvalue, multiplicity, eigenvector" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "ev1 = Matrix([-1, 1])\n", "ev1" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "M * ev1" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "SymPy can also find roots (including their multiplicity), solve systems of linear and non-linear equations and differential equations.\n", "\n", "For more information, see the [SymPy documentation](https://docs.sympy.org/latest/index.html), and to get an overview of the possibilities, take a look at the [SymPy Modules Reference](https://docs.sympy.org/latest/modules/index.html)." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Harmonic oscillator in SymPy and SciPy\n", "As an easy example we will take a look at the differential equation that describes a harmonic oscillator and solve it numerically with SciPy and symbolically with SymPy: \n", "$$\\ddot x + \\omega_0^2 x = 0 \\quad \\text{(differential equation of 2nd order)}$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### SciPy\n", "We start with the numerical solution in SciPy. As with the free-fall example, we rewrite the DE of 2nd order as a system of two coupled DEs of 1st order." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "from scipy.integrate import odeint\n", "\n", "def F(y, t):\n", " # define system of ODEs\n", " omega0 = 1\n", " return [y[1], -omega0**2 * y[0]]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "What `F` actually defines can be written in this form:\n", "$$\\frac{\\mathrm d}{\\mathrm dt} \\left( \\begin{matrix}x \\\\ \\dot x\\end{matrix} \\right) \n", " = \\left( \\begin{matrix} \\dot x \\\\ -\\omega_0^2 x \\end{matrix} \\right)$$" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# array of time values to study\n", "t_min = 0\n", "t_max = 10\n", "dt = 0.1\n", "ts = np.arange(t_min, t_max+dt, dt)\n", "\n", "# initial conditions (x0, v0)\n", "y0 = (1.0, 0.0)\n", "\n", "# solve...\n", "y = odeint(F, y0, ts)\n", "# y[:,0] corresponds to position (=height)\n", "# y[:,1] corresponds to velocity" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### SymPy\n", "Before plotting the result, we also determine the analytical solution with SymPy, using Sympy's `dsolve` to solve the differential equation." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "from sympy import symbols, Function, dsolve\n", "x = symbols('x', cls=Function) # create an undefined function\n", "t, omega0 = symbols('t omega_0') # symbols for time and frequency (note the LaTeX notation as argument)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "omega0" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# define the DE\n", "DE = x(t).diff(t).diff(t) + omega0**2 * x(t)\n", "DE" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# ...and solve it\n", "sol_x = dsolve(DE, x(t))\n", "# 2nd order ODE => two independent solutions\n", "# ODE = ordinary differential equation, i.e. dependent on single independent variable\n", "sol_x" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now we use the right-hand side (`rhs`) of the solution (`sol_x`) and our initial conditions from above ($x(0) = x_0$, $\\dot x(0) = v_0$) to obtain relations from which we can determine $C_1$ and $C_2$." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "from sympy import Eq, solve\n", "x0, v0 = symbols(\"x_0 v_0\")\n", "initial1 = Eq(sol_x.rhs.subs(t, 0), x0)\n", "initial1" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "initial2 = Eq(sol_x.rhs.diff(t).subs(t, 0), v0)\n", "initial2" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "And we solve these two equations to get $C_1$ and $C_2$" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "Cs = solve([initial1, initial2])\n", "Cs " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "which we substitute in the solution for $x(t)$:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "sol_xs = sol_x.subs(Cs[0]) # there is only one solution, so we pick that with [0]\n", "sol_xs = simplify(sol_xs)\n", "sol_xs" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now is this what we would expect? Let us insert the initial conditions from above:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "xt = simplify(sol_xs.subs({x0: 1, v0: 0, omega0: 1}))\n", "xt" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We also define the velocity:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "vt = Eq(xt.lhs.diff(t), xt.rhs.diff(t)) # xt.diff(t) doesn't evaluate\n", "vt" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Plotting the results\n", "SymPy has its own plotting function:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "from sympy import plot\n", "p1 = plot(xt.rhs, show=False, xlim = [t_min, t_max], legend = True)\n", "p2 = plot(vt.rhs, show=False)\n", "p1.append(p2[0])\n", "p1[0].label = xt.lhs\n", "p1[1].label = \"v(t)\"\n", "p1[0].line_color = 'b'\n", "p1[1].line_color = 'r'\n", "p1.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "To plot the numerical result, we use `matplotlib` as before:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "%matplotlib inline\n", "import matplotlib.pyplot as plt\n", "plt.plot(ts, y[:, 0], \"b.\", linewidth=2, label = \"position\")\n", "plt.plot(ts, y[:, 1], \"r.\", linewidth=2, label = \"velocity\");" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Try it out\n", "* Compute the [Padé approximant](https://en.wikipedia.org/wiki/Pad%C3%A9_approximant) of the sine function to third order and compare its behavior to the Taylor series of same order.\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "@webio": { "lastCommId": null, "lastKernelId": null }, "kernelspec": { "display_name": "Python 3 (ipykernel)", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.8.10" } }, "nbformat": 4, "nbformat_minor": 2 }