book/12-symbolic-computation.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Symbolic computation\n",
    "\n",
    "## SymPy\n",
    "\n",
    "In this section, we introduce some basic functionality of the SymPy (SYMbolic Python) library. In contrast to numerical computation (involving numbers), in symbolic calculation we are processing and transforming generic variables.\n",
    "\n",
    "The SymPy home page is <https://www.sympy.org/>, and provides the full (and up-to-date) documentation for this library.\n",
    "\n",
    "Symbolic calculation is very slow compared to floating point operation (see for example [symbolic calculation for decimals](13 Numeric Computation, Symbolic Calculation)), and thus generally not for direct simulation. However, it is a powerful tool to support the preparation of code and symbolic work. Occasionally, we use symbolic operations in simulations to work out the most efficient numerical code, before that is executed.\n",
    "\n",
    "### Output\n",
    "\n",
    "Before we start using sympy, we'll call `init_printing`. This tells sympy to display expressions in a nicer format."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import sympy\n",
    "sympy.init_printing(use_latex='mathjax')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Symbols\n",
    "\n",
    "Before we can carry out any symbolic operations, we need to create symbolic variables using SymPy’s `Symbol` function:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "sympy.core.symbol.Symbol"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from sympy import Symbol\n",
    "x = Symbol('x')\n",
    "type(x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle x$"
      ],
      "text/plain": [
       "x"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "y = Symbol('y')\n",
    "2 * x - x"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle 2 x + 11 y$"
      ],
      "text/plain": [
       "2⋅x + 11⋅y"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x + y + x + 10*y"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle x + 10$"
      ],
      "text/plain": [
       "x + 10"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "y + x - y + 10"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can abbreviate the creation of multiple symbolic variables using the `symbols` function. For example, to create the symbolic variables `x`, `y` and `z`, we can use"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle 2 y + 3 z$"
      ],
      "text/plain": [
       "2⋅y + 3⋅z"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import sympy\n",
    "x, y, z = sympy.symbols('x,y,z')\n",
    "x + 2*y + 3*z - x"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Once we have completed our term manipulation, we sometimes like to insert numbers for variables. This can be done using the `subs` method."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle x + 2 y$"
      ],
      "text/plain": [
       "x + 2⋅y"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from sympy import symbols\n",
    "x, y = symbols('x,y')\n",
    "x + 2*y"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle x + 2 y$"
      ],
      "text/plain": [
       "x + 2⋅y"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x + 2*y.subs(x, 10)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle 2 y + 10$"
      ],
      "text/plain": [
       "2⋅y + 10"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "(x + 2*y).subs(x, 10)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle 16$"
      ],
      "text/plain": [
       "16"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "(x + 2*y).subs(x, 10).subs(y, 3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle 16$"
      ],
      "text/plain": [
       "16"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "(x + 2*y).subs({x:10, y:3})"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can also substitute a symbolic variable for another one such as in this example where `y` is replaced with `x` before we substitute `x` with the number `2`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle 3 x + y^{2}$"
      ],
      "text/plain": [
       "       2\n",
       "3⋅x + y "
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "myterm = 3*x + y**2\n",
    "myterm"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle y^{2} + 3 y$"
      ],
      "text/plain": [
       " 2      \n",
       "y  + 3⋅y"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "myterm.subs(x, y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle 10$"
      ],
      "text/plain": [
       "10"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "myterm.subs(x, y).subs(y, 2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "From this point onward, some of the code fragments and examples we present will assume that the required symbols have already been defined. If you try an example and SymPy gives a message like `NameError: name ’x’ is not defined` it is probably because you need to define the symbol using one of the methods above."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### isympy\n",
    "\n",
    "The `isympy` executable is a wrapper around ipython which creates the symbolic (real) variables `x`, `y` and `z`, the symbolic integer variables `k`, `m` and `n` and the symbolic function variables `f`, `g` and `h`, and imports all objects from the SymPy toplevel.\n",
    "\n",
    "This is convenient to figure out new features or experimenting interactively\n",
    "\n",
    "    $> isympy\n",
    "    Python 2.6.5 console for SymPy 0.6.7\n",
    "\n",
    "    These commands were executed:\n",
    "    >>> from __future__ import division\n",
    "    >>> from sympy import *\n",
    "    >>> x, y, z = symbols('xyz')\n",
    "    >>> k, m, n = symbols('kmn', integer=True)\n",
    "    >>> f, g, h = map(Function, 'fgh')\n",
    "\n",
    "    Documentation can be found at https://www.sympy.org/\n",
    "\n",
    "    In [1]: "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Numeric types\n",
    "\n",
    "SymPy has the numeric types `Rational` and `RealNumber`. The Rational class represents a rational number as a pair of two integers: the numerator and the denominator, so `Rational(1,2)` represents `1/2`, `Rational(5,2)` represents `5/2` and so on."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sympy import Rational"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\frac{1}{10}$"
      ],
      "text/plain": [
       "1/10"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "a = Rational(1, 10)\n",
    "a"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\frac{45}{67}$"
      ],
      "text/plain": [
       "45\n",
       "──\n",
       "67"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "b = Rational(45, 67)\n",
    "b"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\frac{9}{134}$"
      ],
      "text/plain": [
       "9/134"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "a * b"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle - \\frac{383}{670}$"
      ],
      "text/plain": [
       "-383 \n",
       "─────\n",
       " 670 "
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "a - b"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\frac{517}{670}$"
      ],
      "text/plain": [
       "517\n",
       "───\n",
       "670"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "a + b"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that the Rational class works with rational expressions *exactly*. This is in contrast to Python’s standard `float` data type which uses floating point representation to *approximate* (rational) numbers.\n",
    "\n",
    "We can convert the `sympy.Rational` type into a Python floating point variable using `float` or the `evalf` method of the Rational object. The `evalf` method can take an argument that specifies how many digits should be computed for the floating point approximation (not all of those may be used by Python’s floating point type of course)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\frac{2}{3}$"
      ],
      "text/plain": [
       "2/3"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "c = Rational(2, 3)\n",
    "c"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle 0.666666666666667$"
      ],
      "text/plain": [
       "0.6666666666666666"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "float(c)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle 0.666666666666667$"
      ],
      "text/plain": [
       "0.666666666666667"
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "c.evalf()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle 0.66666666666666666666666666666666666666666666666667$"
      ],
      "text/plain": [
       "0.66666666666666666666666666666666666666666666666667"
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "c.evalf(50)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Differentiation and Integration\n",
    "\n",
    "SymPy is capable of carrying out differentiation and integration of many functions:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\cos{\\left(x \\right)}$"
      ],
      "text/plain": [
       "cos(x)"
      ]
     },
     "execution_count": 25,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from sympy import Symbol, exp, sin, sqrt, diff\n",
    "x = Symbol('x')\n",
    "y = Symbol('y')\n",
    "diff(sin(x), x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle 0$"
      ],
      "text/plain": [
       "0"
      ]
     },
     "execution_count": 26,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "diff(sin(x), y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle 9 x^{8} + 20 x + 3$"
      ],
      "text/plain": [
       "   8           \n",
       "9⋅x  + 20⋅x + 3"
      ]
     },
     "execution_count": 27,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "diff(10 + 3*x + 4*y + 10*x**2 + x**9, x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle 4$"
      ],
      "text/plain": [
       "4"
      ]
     },
     "execution_count": 28,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "diff(10 + 3*x + 4*y + 10*x**2 + x**9, y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle 32$"
      ],
      "text/plain": [
       "32"
      ]
     },
     "execution_count": 29,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "diff(10 + 3*x + 4*y + 10*x**2 + x**9, x).subs(x,1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle 263.66015625$"
      ],
      "text/plain": [
       "263.660156250000"
      ]
     },
     "execution_count": 30,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "diff(10 + 3*x + 4*y + 10*x**2 + x**9, x).subs(x,1.5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle e^{x}$"
      ],
      "text/plain": [
       " x\n",
       "ℯ "
      ]
     },
     "execution_count": 31,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "diff(exp(x), x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle - x e^{- \\frac{x^{2}}{2}}$"
      ],
      "text/plain": [
       "      2 \n",
       "    -x  \n",
       "    ────\n",
       "     2  \n",
       "-x⋅ℯ    "
      ]
     },
     "execution_count": 32,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "diff(exp(-x ** 2 / 2), x)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The SymPy `diff()` function takes a minimum of two arguments: the function to be differentiated and the variable with respect to which the differentiation is performed. Higher derivatives may be calculated by specifying additional variables, or by adding an optional integer argument:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle 12 x^{3}$"
      ],
      "text/plain": [
       "    3\n",
       "12⋅x "
      ]
     },
     "execution_count": 33,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "diff(3*x**4, x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle 72 x$"
      ],
      "text/plain": [
       "72⋅x"
      ]
     },
     "execution_count": 34,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "diff(3*x**4, x, x, x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle 72 x$"
      ],
      "text/plain": [
       "72⋅x"
      ]
     },
     "execution_count": 35,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "diff(3*x**4, x, 3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle 1512 x^{2} y^{5}$"
      ],
      "text/plain": [
       "      2  5\n",
       "1512⋅x ⋅y "
      ]
     },
     "execution_count": 36,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "diff(3*x**4*y**7, x, 2, y, 2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle 1512 x^{2} y^{5}$"
      ],
      "text/plain": [
       "      2  5\n",
       "1512⋅x ⋅y "
      ]
     },
     "execution_count": 37,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "diff(diff(3*x**4*y**7, x, x), y, y)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "At times, SymPy may return a result in an unfamiliar form. If, for example, you wish to use SymPy to check that you differentiated something correctly, a technique that might be of use is to subtract the SymPy result from your result, and check that the answer is zero.\n",
    "\n",
    "Taking the simple example of a multiquadric radial basis function, $\\phi(r)=\\sqrt{r^2+\\sigma^2}$ with $r=\\sqrt{x^2+y^2}$ and *σ* a constant, we can verify that the first derivative in *x* is $\\partial\\phi/\\partial x=x/\\sqrt{r^2+\\sigma^2}$.\n",
    "\n",
    "In this example, we first ask SymPy to print the derivative. See that it is printed in a different form to our trial derivative, but the subtraction verifies that they are identical:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "x/sqrt(x**2 + y**2 + σ**2)\n"
     ]
    }
   ],
   "source": [
    "r = sqrt(x**2 + y**2)\n",
    "sigma = Symbol('σ')\n",
    "def phi(x,y,sigma):\n",
    "    return sqrt(x**2 + y**2 + sigma**2)\n",
    "\n",
    "mydfdx= x / sqrt(r**2 + sigma**2)\n",
    "print(diff(phi(x, y, sigma), x))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0\n"
     ]
    }
   ],
   "source": [
    "print(mydfdx - diff(phi(x, y, sigma), x))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here it is trivial to tell that the expressions are identical without SymPy’s help, but in more complicated examples there may be many more terms and it would become increasingly difficult, time consuming and error-prone to attempt to rearrange our trial derivative and SymPy’s answer into the same form. It is in such cases that this subtraction technique is of most use.\n",
    "\n",
    "Integration uses a similar syntax. For the indefinite case, specify the function and the variable with respect to which the integration is performed:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\frac{x^{3}}{3}$"
      ],
      "text/plain": [
       " 3\n",
       "x \n",
       "──\n",
       "3 "
      ]
     },
     "execution_count": 40,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from sympy import integrate\n",
    "integrate(x**2, x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle x^{2} y$"
      ],
      "text/plain": [
       " 2  \n",
       "x ⋅y"
      ]
     },
     "execution_count": 41,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "integrate(x**2, y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle y \\sin{\\left(x \\right)}$"
      ],
      "text/plain": [
       "y⋅sin(x)"
      ]
     },
     "execution_count": 42,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "integrate(sin(x), y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle - \\cos{\\left(x \\right)}$"
      ],
      "text/plain": [
       "-cos(x)"
      ]
     },
     "execution_count": 43,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "integrate(sin(x), x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle e^{- \\frac{x^{2}}{2}}$"
      ],
      "text/plain": [
       "   2 \n",
       " -x  \n",
       " ────\n",
       "  2  \n",
       "ℯ    "
      ]
     },
     "execution_count": 44,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "integrate(-x*exp(-x**2/2), x)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can calculate definite integrals by providing `integrate()` with a tuple containing the variable of interest, the lower and the upper bounds. If several variables are specified, multiple integration is performed. When SymPy returns a result in the `Rational` class, it is possible to evaluate it to a floating-point representation at any desired precision (see [numeric types](#Numeric-types))."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle 1$"
      ],
      "text/plain": [
       "1"
      ]
     },
     "execution_count": 45,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "integrate(x*2, (x, 0, 1))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\frac{x^{3}}{3}$"
      ],
      "text/plain": [
       " 3\n",
       "x \n",
       "──\n",
       "3 "
      ]
     },
     "execution_count": 46,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "integrate(x**2, x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\frac{x^{4}}{12}$"
      ],
      "text/plain": [
       " 4\n",
       "x \n",
       "──\n",
       "12"
      ]
     },
     "execution_count": 47,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "integrate(x**2, x, x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 48,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\frac{x^{4} y}{12}$"
      ],
      "text/plain": [
       " 4  \n",
       "x ⋅y\n",
       "────\n",
       " 12 "
      ]
     },
     "execution_count": 48,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "integrate(x**2, x, x, y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\frac{8}{3}$"
      ],
      "text/plain": [
       "8/3"
      ]
     },
     "execution_count": 49,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "integrate(x**2, (x, 0, 2))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 50,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\frac{16}{3}$"
      ],
      "text/plain": [
       "16/3"
      ]
     },
     "execution_count": 50,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "integrate(x**2, (x, 0, 2), (x, 0, 2), (y, 0, 1))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 51,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle 2.66666666666667$"
      ],
      "text/plain": [
       "2.6666666666666665"
      ]
     },
     "execution_count": 51,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "float(integrate(x**2, (x, 0, 2)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 52,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "sympy.core.numbers.Rational"
      ]
     },
     "execution_count": 52,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "type(integrate(x**2, (x, 0, 2)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 53,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle 2.66666666666667$"
      ],
      "text/plain": [
       "2.66666666666667"
      ]
     },
     "execution_count": 53,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "result_rational=integrate(x**2, (x, 0, 2))\n",
    "result_rational.evalf()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 54,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle 2.6666666666666666666666666666666666666666666666667$"
      ],
      "text/plain": [
       "2.6666666666666666666666666666666666666666666666667"
      ]
     },
     "execution_count": 54,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "result_rational.evalf(50)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Ordinary differential equations\n",
    "\n",
    "SymPy has inbuilt support for solving several kinds of ordinary differential equation via its `dsolve` command. We need to set up the ODE and pass it as the first argument, `eq`. The second argument is the function `f(x)` to solve for. An optional third argument, `hint`, influences the method that `dsolve` uses: some methods are better-suited to certain classes of ODEs, or will express the solution more simply, than others.\n",
    "\n",
    "To set up the ODE solver, we need a way to refer to the unknown function for which we are solving, as well as its derivatives. The `Function` and `Derivative` classes facilitate this:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 55,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle y{\\left(x \\right)} = C_{1} e^{- 5 x}$"
      ],
      "text/plain": [
       "           -5⋅x\n",
       "y(x) = C₁⋅ℯ    "
      ]
     },
     "execution_count": 55,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from sympy import Symbol, dsolve, Function, Derivative, Eq\n",
    "y = Function(\"y\")\n",
    "x = Symbol('x')\n",
    "y_ = Derivative(y(x), x)\n",
    "dsolve(y_ + 5*y(x), y(x))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note how `dsolve` has introduced a constant of integration, `C1`. It will introduce as many constants as are required, and they will all be named `Cn`, where `n` is an integer. Note also that the first argument to `dsolve` is taken to be equal to zero unless we use the `Eq()` function to specify otherwise:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 56,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle y{\\left(x \\right)} = C_{1} e^{- 5 x}$"
      ],
      "text/plain": [
       "           -5⋅x\n",
       "y(x) = C₁⋅ℯ    "
      ]
     },
     "execution_count": 56,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "dsolve(y_ + 5*y(x), y(x))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 57,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle y{\\left(x \\right)} = C_{1} e^{- 5 x}$"
      ],
      "text/plain": [
       "           -5⋅x\n",
       "y(x) = C₁⋅ℯ    "
      ]
     },
     "execution_count": 57,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "dsolve(Eq(y_ + 5*y(x), 0), y(x))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 58,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle y{\\left(x \\right)} = C_{1} e^{- 5 x} + \\frac{12}{5}$"
      ],
      "text/plain": [
       "           -5⋅x   12\n",
       "y(x) = C₁⋅ℯ     + ──\n",
       "                  5 "
      ]
     },
     "execution_count": 58,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "dsolve(Eq(y_ + 5*y(x), 12), y(x))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The results from `dsolve` are an instance of the `Equality` class. This has consequences when we wish to numerically evaluate the function and use the result elsewhere (*e.g.* if we wanted to plot *y*(*x*) against *x*), because even after using `subs()` and `evalf()`, we still have an `Equality`, not any sort of scalar. The way to evaluate the function to a number is via the `rhs` attribute of the `Equality`.\n",
    "\n",
    "Note that, here, we use `z` to store the `Equality` returned by `dsolve`, even though it is an expression for a function called `y(x)`, to emphasise the distinction between the `Equality` itself and the data that it contains."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 59,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle y{\\left(x \\right)} = C_{1} e^{- 5 x}$"
      ],
      "text/plain": [
       "           -5⋅x\n",
       "y(x) = C₁⋅ℯ    "
      ]
     },
     "execution_count": 59,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "z = dsolve(y_ + 5*y(x), y(x))\n",
    "z"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 60,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "sympy.core.relational.Equality"
      ]
     },
     "execution_count": 60,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "type(z)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 61,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle C_{1} e^{- 5 x}$"
      ],
      "text/plain": [
       "    -5⋅x\n",
       "C₁⋅ℯ    "
      ]
     },
     "execution_count": 61,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "z.rhs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 62,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle y{\\left(3 \\right)} = \\frac{2}{e^{15}}$"
      ],
      "text/plain": [
       "          -15\n",
       "y(3) = 2⋅ℯ   "
      ]
     },
     "execution_count": 62,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "C1=Symbol('C1')\n",
    "y3 = z.subs({C1:2, x:3})\n",
    "y3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 63,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle y{\\left(3 \\right)} = 6.11804641 \\cdot 10^{-7}$"
      ],
      "text/plain": [
       "y(3) = 6.11804641e-7"
      ]
     },
     "execution_count": 63,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "y3.evalf(10)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 64,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\frac{2}{e^{15}}$"
      ],
      "text/plain": [
       "   -15\n",
       "2⋅ℯ   "
      ]
     },
     "execution_count": 64,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "y3.rhs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 65,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle 6.11804641 \\cdot 10^{-7}$"
      ],
      "text/plain": [
       "6.118046410e-7"
      ]
     },
     "execution_count": 65,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "y3.rhs.evalf(10)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 66,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle 4.122307245 \\cdot 10^{-9}$"
      ],
      "text/plain": [
       "4.122307245e-9"
      ]
     },
     "execution_count": 66,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "z.rhs.subs({C1:2, x:4}).evalf(10)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 67,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle 2.777588773 \\cdot 10^{-11}$"
      ],
      "text/plain": [
       "2.777588773e-11"
      ]
     },
     "execution_count": 67,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "z.rhs.subs({C1:2, x:5}).evalf(10)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 68,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "sympy.core.numbers.Float"
      ]
     },
     "execution_count": 68,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "type(z.rhs.subs({C1:2, x:5}).evalf(10))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "At times, `dsolve` may return too general a solution. One example is when there is a possibility that some coefficients may be complex. If we know that, for example, they are always real and positive, we can provide `dsolve` this information to avoid the solution becoming unnecessarily complicated:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 69,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle f{\\left(x \\right)} = C_{1} e^{- i a^{2} x} + C_{2} e^{i a^{2} x}$"
      ],
      "text/plain": [
       "               2            2  \n",
       "           -ⅈ⋅a ⋅x       ⅈ⋅a ⋅x\n",
       "f(x) = C₁⋅ℯ        + C₂⋅ℯ      "
      ]
     },
     "execution_count": 69,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from sympy import *\n",
    "a, x = symbols('a,x')\n",
    "f = Function('f')\n",
    "dsolve(Derivative(f(x), x, 2) + a**4*f(x), f(x))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 70,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle f{\\left(x \\right)} = C_{1} \\sin{\\left(a^{2} x \\right)} + C_{2} \\cos{\\left(a^{2} x \\right)}$"
      ],
      "text/plain": [
       "             ⎛ 2  ⎞         ⎛ 2  ⎞\n",
       "f(x) = C₁⋅sin⎝a ⋅x⎠ + C₂⋅cos⎝a ⋅x⎠"
      ]
     },
     "execution_count": 70,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "a = Symbol('a',real=True,positive=True)\n",
    "dsolve(Derivative(f(x), x, 2)+a**4*f(x), f(x))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Series expansions and plotting\n",
    "\n",
    "It is possible to expand many SymPy expressions as Taylor series. The `series` method makes this straightforward. At minimum, we must specify the expression and the variable in which to expand it. Optionally, we can also specify the point around which to expand, the maximum term number, and the direction of the expansion (try `help(Basic.series)` for more information)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 71,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle x - \\frac{x^{3}}{6} + \\frac{x^{5}}{120} + O\\left(x^{6}\\right)$"
      ],
      "text/plain": [
       "     3     5        \n",
       "    x     x     ⎛ 6⎞\n",
       "x - ── + ─── + O⎝x ⎠\n",
       "    6    120        "
      ]
     },
     "execution_count": 71,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from sympy import *\n",
    "x = Symbol('x')\n",
    "sin(x).series(x, 0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 72,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle x - \\frac{x^{3}}{6} + \\frac{x^{5}}{120} + O\\left(x^{6}\\right)$"
      ],
      "text/plain": [
       "     3     5        \n",
       "    x     x     ⎛ 6⎞\n",
       "x - ── + ─── + O⎝x ⎠\n",
       "    6    120        "
      ]
     },
     "execution_count": 72,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "series(sin(x), x, 0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 73,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle 1.11729533119247 - 0.438791280945186 \\left(x - 0.5\\right)^{2} + 0.0799042564340338 \\left(x - 0.5\\right)^{3} + 0.0365659400787655 \\left(x - 0.5\\right)^{4} - 0.00399521282170169 \\left(x - 0.5\\right)^{5} - 0.00121886466929218 \\left(x - 0.5\\right)^{6} + 9.51241148024212 \\cdot 10^{-5} \\left(x - 0.5\\right)^{7} + 2.17654405230747 \\cdot 10^{-5} \\left(x - 0.5\\right)^{8} - 1.32116826114474 \\cdot 10^{-6} \\left(x - 0.5\\right)^{9} - 0.479425538604203 x + O\\left(\\left(x - \\frac{1}{2}\\right)^{10}; x\\rightarrow \\frac{1}{2}\\right)$"
      ],
      "text/plain": [
       "                                              2                               3                               4      \n",
       "1.11729533119247 - 0.438791280945186⋅(x - 0.5)  + 0.0799042564340338⋅(x - 0.5)  + 0.0365659400787655⋅(x - 0.5)  - 0.0\n",
       "\n",
       "                          5                                6                                7                        \n",
       "0399521282170169⋅(x - 0.5)  - 0.00121886466929218⋅(x - 0.5)  + 9.51241148024212e-5⋅(x - 0.5)  + 2.17654405230747e-5⋅(\n",
       "\n",
       "        8                                9                          ⎛         10         ⎞\n",
       "x - 0.5)  - 1.32116826114474e-6⋅(x - 0.5)  - 0.479425538604203⋅x + O⎝(x - 1/2)  ; x → 1/2⎠"
      ]
     },
     "execution_count": 73,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# NBVAL_IGNORE_OUTPUT\n",
    "cos(x).series(x, 0.5, 10)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In some cases, especially for numerical evaluation and plotting the results, it is necessary to remove the trailing `O(n)` term:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 74,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle - 0.479425538604203 x - 1.32116826114474 \\cdot 10^{-6} \\left(x - 0.5\\right)^{9} + 2.17654405230747 \\cdot 10^{-5} \\left(x - 0.5\\right)^{8} + 9.51241148024212 \\cdot 10^{-5} \\left(x - 0.5\\right)^{7} - 0.00121886466929218 \\left(x - 0.5\\right)^{6} - 0.00399521282170169 \\left(x - 0.5\\right)^{5} + 0.0365659400787655 \\left(x - 0.5\\right)^{4} + 0.0799042564340338 \\left(x - 0.5\\right)^{3} - 0.438791280945186 \\left(x - 0.5\\right)^{2} + 1.11729533119247$"
      ],
      "text/plain": [
       "                                                    9                                8                               \n",
       "-0.479425538604203⋅x - 1.32116826114474e-6⋅(x - 0.5)  + 2.17654405230747e-5⋅(x - 0.5)  + 9.51241148024212e-5⋅(x - 0.5\n",
       "\n",
       " 7                                6                                5                               4                 \n",
       ")  - 0.00121886466929218⋅(x - 0.5)  - 0.00399521282170169⋅(x - 0.5)  + 0.0365659400787655⋅(x - 0.5)  + 0.079904256434\n",
       "\n",
       "              3                              2                   \n",
       "0338⋅(x - 0.5)  - 0.438791280945186⋅(x - 0.5)  + 1.11729533119247"
      ]
     },
     "execution_count": 74,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# NBVAL_IGNORE_OUTPUT\n",
    "cos(x).series(x, 0.5, 10).removeO()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "SymPy provides two inbuilt plotting functions, `Plot()` from the `sympy.plotting` module, and `plot` from `sympy.mpmath.visualization`. At the time of writing, these functions lack the ability to add a key to the plot, which means they are unsuitable for most of our needs. Should you wish to use them nevertheless, their `help()` text is useful.\n",
    "\n",
    "For most of our purposes, Matplotlib should be the plotting tool of choice. The details are in chapter \\[cha:visualisingdata\\]. Here we furnish just one example of how to plot the results of a SymPy computation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 75,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 76,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<Figure size 640x480 with 0 Axes>"
      ]
     },
     "execution_count": 76,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "text/plain": [
       "<Figure size 640x480 with 0 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "from sympy import sin,series,Symbol\n",
    "import pylab\n",
    "x = Symbol('x')\n",
    "s10 = sin(x).series(x,0,10).removeO()\n",
    "s20 = sin(x).series(x,0,20).removeO()\n",
    "s = sin(x)\n",
    "xx = []\n",
    "y10 = []\n",
    "y20 = []\n",
    "y = []\n",
    "for i in range(1000):\n",
    "  xx.append(i / 100.0)\n",
    "  y10.append(float(s10.subs({x:i/100.0})))\n",
    "  y20.append(float(s20.subs({x:i/100.0})))\n",
    "  y.append(float(s.subs({x:i/100.0})))\n",
    "\n",
    "pylab.figure()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 77,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x7fea785d47f0>"
      ]
     },
     "execution_count": 77,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAj4AAAG2CAYAAAB/OYyEAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8o6BhiAAAACXBIWXMAAA9hAAAPYQGoP6dpAABjuElEQVR4nO3dd1zV1ePH8de9l42IAiqigrhyj7Rh5cqV2rBhVpaaLX9ZZlZWWmnT1teWpVmullnZ0ErTLDVzYu69xYGoICDIBe79/P64QpFWgBc+3Hvfz8eDx71c7ngDyX13Puecj8UwDAMRERERH2A1O4CIiIhIWVHxEREREZ+h4iMiIiI+Q8VHREREfIaKj4iIiPgMFR8RERHxGSo+IiIi4jNUfERERMRnqPiIiIiIz1DxEREREZ/hscVn7NixWCwWhg0bZnYUERER8RAeWXxWr17NpEmTaN68udlRRERExIN4XPE5deoU/fr144MPPqBy5cpmxxEREREP4md2gOIaMmQIvXr1okuXLrzwwgv/el+73Y7dbi/43Ol0kpKSQmRkJBaLpbSjioiIiBsYhkFGRgYxMTFYrec3ZuNRxefzzz/njz/+YPXq1UW6/9ixY3n22WdLOZWIiIiUhcTERGrWrHlez+ExxScxMZGHHnqI+fPnExQUVKTHPPnkkwwfPrzg87S0NGJjY0lMTKRixYqlFVVEfFx2roM2L/wMwJIRnYgIDTA5kYhnS09Pp1atWoSFhZ33c3lM8VmzZg3Jycm0bt264DaHw8GSJUsYP348drsdm81W6DGBgYEEBgae9VwVK1ZU8RGRUnPieCbWwBCC/W3ERevQuoi7uOPfkscUn86dO7Nx48ZCt9155500bNiQxx9//KzSIyJilsMnTwNQo3KwSo9IOeMxxScsLIymTZsWui00NJTIyMizbhcRMdOhVFfxiakUbHISEfk7j1vOLiJS3h3MH/FR8REpdzxmxOdcFi1aZHYEEZGz5B/qqllZxcfbOBwOcnNzzY7hdfz9/ctsyopHFx8RkfLoz0NdRVuBKuWfYRgkJSVx8uRJs6N4rUqVKhEdHV3q8+JUfERE3OxQwaGuEJOTiLvkl56qVasSEhKiSetuZBgGWVlZJCcnA1C9evVSfT0VHxERN3I6DY6k/bmqSzyfw+EoKD2RkZFmx/FKwcGufyvJyclUrVq1VA97aXKziIgbHTtlJ9dhYLNaqBZ29j5i4nny5/SEhGgErzTl/3xLew6Vio+IiBsdPDO/J7piEH42/Yn1Jjq8VbrK6uerf5UiIm50WEvZRco1FR8RETfKn9isFV3iibZv3050dDQZGRlufd5HH32UoUOHuvU5S0rFR0TEjfKXsmtis5QniYmJ3HXXXcTExBAQEEBcXBwPPfQQJ06cKHS/UaNGMWTIkIKTgWZnZzNw4ECaNWuGn58fvXv3PufzL168mNatWxMUFESdOnWYOHFioa+PGDGCqVOnsnfv3lL5/opDxUdExI0Oaym7lDN79uyhTZs27NixgxkzZrBr1y4mTpzIwoULadu2LSkpKQAcPHiQ2bNnc+eddxY81uFwEBwczNChQ+nSpcs5n3/v3r307NmTdu3asXbtWkaOHMnQoUOZNWtWwX2qVq1Kt27dzipEZlDxERFxIx3qkvJmyJAhBAQEMH/+fDp06EBsbCw9evTg559/5tChQ4waNQqAL774ghYtWlCzZs2Cx4aGhjJhwgTuueceoqOjz/n8EydOJDY2ljfffJNGjRpx9913M2jQIF5//fVC97v22muZMWNG6X2jRaTiIyLiRvmHunS6Cu9mGAZZOXmmfBiGUeScKSkp/PTTT9x///0Fe+Xki46Opl+/fsycORPDMFiyZAlt2rQp9s9i+fLldOvWrdBt3bt3JyEhodDS9IsvvpjExET2799f7NdwJ21gKCLiJmmnc8mw5wE6M7u3O53roPEzP5ny2lue605IQNHevnfu3IlhGDRq1OicX2/UqBGpqakcO3aMffv20bp162LnSUpKolq1aoVuq1atGnl5eRw/frxgJ+YaNWoAsG/fPuLi4or9Ou6iER8RETdJTMkCIKpCQJHfmETMlD96ZLFYOH36NEFBJTtE+/c9eP76vPnyR5yysrJK9Bruon+ZIiJukl98albWxGZvF+xvY8tz3U177aKqV68eFouFLVu2nHNF1rZt26hcuTJRUVFERUWRmppa7DzR0dEkJSUVui05ORk/P79Cp/jIn0RdpUqVYr+GO6n4iIi4SWKqq/jUilDx8XYWi8UjRvUiIyPp2rUr7733Hg8//HCheT5JSUl8+umn9O/fH4vFQqtWrdiyZUuxX6Nt27bMmTOn0G3z58+nTZs2+Pv7F9y2adMm/P39adKkScm/ITfQoS4RETdJTHFNbI6N0PweKT/Gjx+P3W6ne/fuLFmyhMTERObNm0fXrl2pUaMGL774IuCakLx8+XIcDkehx2/ZsoV169aRkpJCWloa69atY926dQVfHzx4MPv372f48OFs3bqVKVOmMHnyZB599NFCz/Pbb7/Rrl27syZZl7XyX1dFRDxEwYiPDnVJOVK/fn0SEhIYM2YMffv25cSJE0RHR9O7d29Gjx5NREQEAD179sTf35+ff/6Z7t3/PIzXs2fPQiuxWrVqBfw5jyc+Pp4ff/yRhx9+mHfffZeYmBjefvttbrzxxkI5ZsyYwbPPPlva3+5/UvEREXGTAyk61CXlU1xcHFOnTv3X+9hsNkaOHMm4ceMKFZ99+/b95/N36NCBP/744x+//sMPP2Cz2bjpppuKnLm0qPiIiLiB02kUnJldIz7iqe69915SU1PJyMgoOG2FO2RmZjJ16lT8/MyvHeYnEBHxAsdO2cnJc2K1QHXt2iweys/Pr2AnZ3e6+eab3f6cJaXJzSIibpC/lD2mUjD+Nv1pFSmv9K9TRMQNCub36DCXSLmm4iMi4gb5S9lraSm7SLmm4iMi4gZayi7iGVR8RETcIFFL2UU8goqPiIgbqPiIeAYVHxGR85ST5+RIejagOT4i5Z2Kj4jIeTp88jSGAUH+VqpUCDQ7jkiJbd++nejoaDIyMtz6vI8++ihDhw5163OWlIqPiMh5yp/YXLNyCBaLxeQ0ImdLTEzkrrvuIiYmhoCAAOLi4njooYc4ceJEofuNGjWKIUOGFOzavGjRIq677jqqV69OaGgoLVu25NNPPz3r+RcvXkzr1q0JCgqiTp06TJw4sdDXR4wYwdSpU9m7d2/pfZNFpOIjInKe8vfwidX8HimH9uzZQ5s2bdixYwczZsxg165dTJw4kYULF9K2bVtSUlIAOHjwILNnz+bOO+8seOyyZcto3rw5s2bNYsOGDQwaNIj+/fszZ86cgvvs3buXnj170q5dO9auXcvIkSMZOnQos2bNKrhP1apV6dat21mFyAw6ZYWIyHkq2MOnsub3SPkzZMgQAgICmD9/PsHBrv9GY2NjadWqFXXr1mXUqFFMmDCBL774ghYtWlCzZs2Cx44cObLQcw0dOpSffvqJb775hmuuuQaAiRMnEhsby5tvvglAo0aNSEhI4PXXXy90hvZrr72Wp59+mldeeaWUv+N/pxEfEZHztP9EJgBxkaEmJ5EyYxiQk2nOh2EUOWZKSgo//fQT999/f0HpyRcdHU2/fv2YOXMmhmGwZMkS2rRp85/PmZaWRkRERMHny5cvp1u3boXu0717dxISEsjNzS247eKLLyYxMZH9+/cXOX9p0IiPiMh52nfCdairdpQOdfmM3Cx4Kcac1x55GAKKVrJ37tyJYRg0atTonF9v1KgRqampHDt2jH379tG6det/fb6vvvqK1atX8/777xfclpSURLVq1Qrdr1q1auTl5XH8+HGqV68OQI0aNQDYt28fcXFxRcpfGjTiIyJyHgzD0IiPeCzjzOiRxWLh9OnTBAUF/eN9Fy1axMCBA/nggw9o0qRJoa/9fVL/X583X/6IU1ZWlluyl5RGfEREzsOxU3aychxYLVBTc3x8h3+Ia+TFrNcuonr16mGxWNiyZQu9e/c+6+vbtm2jcuXKREVFERUVRWpq6jmfZ/HixVxzzTWMGzeO/v37F/padHQ0SUlJhW5LTk7Gz8+PyMjIgtvyJ1FXqVKlyPlLg8eM+EyYMIHmzZtTsWJFKlasSNu2bZk7d67ZsUTEx+0/c5grplIwgX42k9NImbFYXIebzPgoxpYJkZGRdO3alffee4/Tp08X+lpSUhKffvopffv2xWKx0KpVK7Zs2XLWcyxatIhevXrx8ssvc++995719bZt27JgwYJCt82fP582bdrg7+9fcNumTZvw9/c/a7SorHlM8alZsyYvv/wyCQkJJCQkcOWVV3LdddexefNms6OJiA/bd9x1mKu2DnNJOTV+/Hjsdjvdu3dnyZIlJCYmMm/ePLp27UqNGjV48cUXAdeE5OXLl+NwOAoem196hg4dyo033khSUhJJSUkFozcAgwcPZv/+/QwfPpytW7cyZcoUJk+ezKOPPloox2+//Ua7du3OmmRd1jym+FxzzTX07NmTBg0a0KBBA1588UUqVKjAihUrzI4mIj4sf8QnLlITm6V8ql+/PgkJCdStW5e+fftSt25d7r33Xjp16sTy5csLVmj17NkTf39/fv7554LHTps2jaysLMaOHUv16tULPm644YaC+8THx/Pjjz+yaNEiWrZsyfPPP8/bb79daCk7wIwZM7jnnnvK5pv+Fx45x8fhcPDll1+SmZlJ27Zt//F+drsdu91e8Hl6enpZxBMRH7LvhEZ8pPyLi4tj6tSp/3ofm83GyJEjGTduHN27dwdcxWfatGn/+fwdOnTgjz/++Mev//DDD9hsNm666aZi5S4NHlV8Nm7cSNu2bcnOzqZChQp88803NG7c+B/vP3bsWJ599tkyTCgivkYjPuJN7r33XlJTU8nIyCg4bYU7ZGZmMnXqVPz8zK8dFsMoxk5IJsvJyeHAgQOcPHmSWbNm8eGHH7J48eJ/LD/nGvGpVasWaWlpVKxYsaxii4iXMgyD5s/OJyM7j/kPt6dBNfe9UUj5kZ2dzd69e4mPj//X5d5yfv7t55yenk54eLhb3r/Nr17FEBAQQL169QBo06YNq1ev5q233iq0kdJfBQYGEhioMyWLSOlIzcolIzsP0Hm6RDyFx0xuPhfDMAqN6IiIlKX8+T3Vw4MI8tdSdhFP4DEjPiNHjqRHjx7UqlWLjIwMPv/8cxYtWsS8efPMjiYiPurPHZs12iPiKTym+Bw9epQ77riDI0eOEB4eTvPmzQv2IRARMUP+xGat6BLxHB5TfCZPnmx2BBGRQv5c0aXiI+IpPHqOj4iImf7cw0eHukQ8hYqPiEgJacRHxPOo+IiIlEBaVi4pmTmAJjeL5xk4cOA5z9b+X3755RcaNmyI0+ks0v03btxIzZo1yczMLPZrlRYVHxGREth9/BQA0RWDCA30mOmSIgC89dZbRToVxd+NGDGCUaNGYbUWrT40a9aMiy++mDfeeKPYr1VaVHxEREpgd7Kr+NSposNc4nnCw8OpVKlSsR6zbNkydu7cSZ8+fYr1uDvvvJMJEyYUOuu7mVR8RERKYM9x19C9io+UZ1999RXNmjUjODiYyMhIunTpQmZm5lmHujp27MjQoUMZMWIEERERREdHM2bMmELP9fnnn9OtW7eC00kYhkGXLl246qqryD/71cmTJ4mNjWXUqFEFj+vevTsnTpxg8eLFpf79FoWKj4hICew5dmbEJ6qCyUnEDIZhkJWbZcpHUU+xeeTIEW699VYGDRrE1q1bWbRoETfccMM/Pn769OmEhoaycuVKXn31VZ577jkWLFhQ8PUlS5bQpk2bgs8tFgvTp09n1apVvP322wAMHjyYatWqFSpNAQEBtGjRgt9++60EP2n304FpEZES2HNMIz6+7HTeaS757BJTXnvlbSsJ8f/vCfVHjhwhLy+PG264gbi4OMA15+afNG/enNGjRwNQv359xo8fz8KFCws2Ct63bx8xMTGFHlOjRg3ef/997rjjDo4ePcqcOXNYu3Yt/v7+Z91v3759xfk2S41GfEREisnhNAqWstetohEfKZ9atGhB586dadasGX369OGDDz4gNTX1H+/fvHnzQp9Xr16d5OTkgs9Pnz59zrPT9+nThxtuuIGxY8fyv//9jwYNGpx1n+DgYLKyss7ju3EfjfiIiBTTwdQschxOAvysxFQKNjuOmCDYL5iVt6007bWLwmazsWDBApYtW8b8+fN55513GDVqFCtXnjv330dpLBZLoWXrUVFR5yxOWVlZrFmzBpvNxs6dO8/53CkpKdStW7dIuUubio+ISDHlH+aKjwzFZrWYnEbMYLFYinS4yWwWi4XLL7+cyy+/nGeeeYa4uDi++eabEj1Xq1at2LJly1m3P/LII1itVubOnUvPnj3p1asXV155ZaH7bNq0iZtuuqlEr+tuOtQlIlJMu49pKbuUfytXruSll14iISGBAwcO8PXXX3Ps2DEaNWpUoufr3r07S5cuLXTbDz/8wJQpU/j000/p2rUrTzzxBAMGDCg0MrRv3z4OHTpEly5dzuv7cRcVHxGRYtJSdvEEFStWZMmSJfTs2ZMGDRrw1FNP8b///Y8ePXqU6Pluv/12tmzZwvbt2wE4duwYd911F2PGjOHCCy8EYPTo0cTExDB48OCCx82YMYNu3boVTLA2m8Uo6ro4L5Cenk54eDhpaWlUrFjR7Dgi4qFumbScFXtSGHdzC264sKbZcaSUZWdns3fvXuLj4885udeXjBgxgrS0NN5///0i3d9ut1O/fn1mzJjB5Zdf/q/3/befszvfvzXiIyJSTH8uZdeKLvEto0aNIi4ursi7MO/fv59Ro0b9Z+kpS5rcLCJSDBnZuSRn2AEd6hLfEx4ezsiRI4t8/wYNGpxzebuZNOIjIlIMe8/M74mqEEjFIP//uLeIlDcqPiIixaAVXSKeTcVHRKQY8uf31FXx8Tk+tBbIFGX181XxEREphp1HXSM+9aqGmZxEykr+jsbl5ZQL3ir/5/v3HaTdTZObRUSKYUdyBgANqmlFl6+w2WxUqlSp4LxVISEhWCzasdtdDMMgKyuL5ORkKlWqhM1mK9XXU/ERESkie56j4OSk9TXi41Oio6MBCp20U9yrUqVKBT/n0qTiIyJSRHuPZ+JwGoQF+lGtYqDZcaQMWSwWqlevTtWqVcnNzTU7jtfx9/cv9ZGefCo+IiJFlD+/p361CjrU4aNsNluZvUFL6dDkZhGRItqZfKb46DCXiMdS8RERKaKdR10Tm+trYrOIx1LxEREpooIRn2oa8RHxVCo+IiJFkJPnZN+Z01XUr6oRHxFPpeIjIlIE+05kkuc0qBDoR/XwILPjiEgJqfiIiBTBnzs2a0WXiCdT8RERKYId+RObdZhLxKOp+IiIFMGu5D/38BERz6XiIyJSBDuT85eya0WXiCdT8RER+Q85eU72akWXiFdQ8RER+Q+7j50i1+E6R1eNSsFmxxGR8+AxxWfs2LFcdNFFhIWFUbVqVXr37s327dvNjiUiPmBbUjoADauHaUWXiIfzmOKzePFihgwZwooVK1iwYAF5eXl069aNzMxMs6OJiJfbdsQ1v6dhdEWTk4jI+fKYs7PPmzev0OdTp06latWqrFmzhvbt25uUSkR8wdakM8WnuiY2i3g6jxnx+bu0tDQAIiIiTE4iIt5u25Ezh7o04iPi8TxmxOevDMNg+PDhXHHFFTRt2vQf72e327Hb7QWfp6enl0U8EfEiJ07ZSc5w/R1pGK0RHxFP55EjPg888AAbNmxgxowZ/3q/sWPHEh4eXvBRq1atMkooIt5i+5nDXHGRIYQGeuT/K4rIX3hc8XnwwQeZPXs2v/76KzVr1vzX+z755JOkpaUVfCQmJpZRShHxFgXzezTaI+IVPOZ/XwzD4MEHH+Sbb75h0aJFxMfH/+djAgMDCQwMLIN0IuKtNL9HxLt4TPEZMmQIn332Gd999x1hYWEkJSUBEB4eTnCwNhQTkdKx7cyITyOt6BLxCh5zqGvChAmkpaXRsWNHqlevXvAxc+ZMs6OJiJfKczgLzsquER8R7+AxIz6GYZgdQUR8zL4TWdjznAT724iNCDE7joi4gceM+IiIlLX8U1VcEB2G1apTVYh4AxUfEZF/sOWwq/hofo+I91DxERH5B5vOFJ+mNcJNTiIi7qLiIyJyDoZhsPmQ69Q4TWNUfES8hYqPiMg5HEnL5kRmDjarhQu0eaGI11DxERE5h01nRnvqV61AkL/N5DQi4i4qPiIi56D5PSLeScVHROQc8uf3NFPxEfEqKj4iIuew6fCZic01tGOziDdR8RER+ZvkjGyOptuxWKBRdRUfEW+i4iMi8jebz8zvqVulAiEBHnNmHxEpAhUfEZG/+XP/Ho32iHgbFR8Rkb/ZdEgrukS8lYqPiMjfbDwz4tNEOzaLeB0VHxGRvzh+ys6hk6exWLSiS8QbqfiIiPzF+sSTANSrUoGwIH9zw4iI26n4iIj8xbozxadFrUqm5hCR0qHiIyLyF/nFp6WKj4hXUvERETnD6TQKDnWp+Ih4JxUfEZEz9p7IJD07j0A/KxdEh5kdR0RKgYqPiMgZ+aM9TWuE42/Tn0cRb6R/2SIiZ2h+j4j3U/ERETlDxUfE+6n4iIgA2bkOth5xnapCxUfEe6n4iIgAW46kk+swiAwNoGblYLPjiEgpUfEREQHWHjgJuDYutFgs5oYRkVKj4iMiAqzZnwJA67jKJicRkdKk4iMiPs8wDFbvSwWgjYqPiFdT8RERn5eYcppjGXb8bRado0vEy6n4iIjPW73PdZirWY1wgvxtJqcRkdKk4iMiPi/hzPyeNrUjTE4iIqVNxUdEfF6C5veI+AwVHxHxaamZOexMPgVoRZeIL1DxERGftma/a7SnTpVQIisEmpxGREqbio+I+LSEM8XnojjN7xHxBSo+IuLTEs6s6GpdW4e5RHyBRxWfJUuWcM011xATE4PFYuHbb781O5KIeLDTOQ7WHzwJwEVa0SXiEzyq+GRmZtKiRQvGjx9vdhQR8QJr9qeS6zCIrhhE7cgQs+OISBnwMztAcfTo0YMePXqYHcNrOR0Ojp1M4kDSNg4f38tpewa5eXZyHXYAggIqEBoUTlhwJWpWa0CdGo3w8/M3ObVIyS3fcxyAy+pG6sSkXs5pOEmzp5GSncJJ+0nsDjsOp4M8Zx5Wi5UQ/xBC/EII8Q+hSnAVKgRUMDuylBKPKj7FZbfbsdvtBZ+np6ebmKb8yMvLZdXmn1m9Yx4H03dxLO8YydYsjvlBtrWIf/y3gp9hEOEwiHT4U80SQY2Q2tSLbk2HC6+nSuWY0v0mRNxg+e4TAFxaN9LkJOIu2XnZbEvZxuYTm9mbtpf96fs5kH6ApKwknIazyM8T6h9KdEg0NcJqUL9SfRpUbkCDyg2ID4/HZtXu3p7Mq4vP2LFjefbZZ82OYTqnw8HC1V+xdPvX7MrexV7/bDJsZ45y+p354M/CE+Q0CHcYBBpW/ACb4fpajsUgx2Jw2mqQZrWQZ7GQ7Gch2c/BVo5B3jE4uBpb4gRq51qJt0TTpGpbrm//f0RWii7rb1vkX52y57H+YBoAbeuo+Hiq03mnSUhKYNnhZaw5uoadqTvJM/L+8f5hAWFUDqxMkF8QflY//Cx+OA0nWXlZZOVlcSrnFKdyT5GZm8nutN3sTtvNkoNL/ny8fxgtq7akdbXWXBx9MU2immC1eNSsEZ9nMQzDMDtESVgsFr755ht69+79j/c514hPrVq1SEtLo2LFimWQ0jxZ2Zl89ctbLDs4l622FFL8Cv/DDHQa1M61Ud0SQbWgWtSOakKDmq2Jr9m4SKM1WdmZ7ErcyL4jm9h3dCOJ6TtJykvmkN9pjv3ttQKcBg1zAmlaoTk3XDaUC+JbufV7FSmJX7cnc+fU1dSsHMzSx680O44UQ0p2CgsPLGTBvgUkHE0g15lb6OsRQRE0i2pG/cr1iQ2LpXZ4bWpUqEHloMr4W//78HxWbhZJWUkkZSaxP30/O1N3sjN1JztSd5CVl1XovlHBUbSv2Z6ONTtyWY3LCLRpL6jSkJ6eTnh4uFvev726+PydO39w5ZHT4WD2kg9ZsOsz1tmOk277s4AEOQ0a5gTRMLQJret0o/2FvQkJCi2VHBt2LGPJxq/ZnrKGbZZjJPn/OZpkNQya2AO4uPIV3Nr5capF1iiVDCL/ZeyPW3l/yR5ublOTV29qYXYc+Q/Zedks2L+AObvnsCppFQ7DUfC16JBqXB7ZjEtCY2lhC6O608BiOMGZBxYbBIVDcCUIrQqRdV3XSyDPmcf21O38cfQPEpISWJm0kszczIKvhwWE0b12d66pcw2tqrbSvDE3UvFBxeevDh/bz/T5Y/g9K4H9AX/eXtnhpKWjGpfFXU2vywcRFlqpzLM5HQ6WbZjLwo2fsCF7KzsC/zzGHuJ0ckleNfq0GU67VleXeTbxbdeOX8qGg2m80bcF17eqaXYc+Qe7Uncxa+csZu+eTXrOn/M0GwdE0s3hT6eTx4k/sc9VdIoqtCpUbQSxl0LcZVDrUvAPKna2XEcuq4+uZnHiYhYeWMjRrKMFX4sNi+XWhrdyXb3rCAsIK/ZzS2E+W3xOnTrFrl27AGjVqhXjxo2jU6dOREREEBsb+5+P97bi88e235j222iW+yUXTEoOdBq0ya1Mx9rX07vDYIICy9cS3dWbF/LtqndZ49jBob+MBDXJ9qNrzHUM6DFKK8Wk1KVn59Ly2fk4DVjxZGeiw4v/pielxzAMlh9ezpRNU1iZtLLg9hj8uOHkSXpmpFMr72/zeALCICwaKlRzjejY/MHq5xr1yU6D0ych/TCcSjr7BQMqQP1u0KQ3NOgBfgFn3+c/OJwOEo4mMGf3HBbsX1BwSCzYL5hr615Lv0b9iA+PL/bziovPFp9FixbRqVOns24fMGAA06ZN+8/He0vxWZTwDZ/+8RqrA9JxnBlKjc2By4JbM6DbGGpWrW1uwCJwOhx8sfAt5u39nLWBWTjPfB+1c6BHZE/uvvo5AgJ0rFxKx89bjnL3RwnER4Xy66MdzY4jZzicDhbsX8CUTVPYmrIVAJsBHbOy6JNxirans12bz4VWgdi2rhGbqo1dozcVqkFRDi1lp8OJXXBkHexfBnt/K1yGQqtAq9vh4nuhYslWp2blZjFn9xxmbJvB7rTdAFiwcFXtq7i3+b3Uq1yvRM/ry3y2+JwvTy8+y9bP5YMVY0gI+nNyXTO7P9fWvoObOw/FavPMJZYbdizjo99e4DfrfrKsrnlJNXMNelbqxv/1fkUjQOJ2z3y3iY+W7+f2S2N5oXczs+P4PMMwWHhgIe+sfYc9aXsACHYa3JiRQf+0DKo7nFDrYrigJ1zQA6IaFK3kFIXTCYf/gC3fwYaZcOrM4SpbAFzYH654GMJLdijUMAxWJq3k0y2fsujgooLbu8Z1ZXCLwTSo3MAN34BvUPEpIU8tPpt3J/DeL4/wu/8JHBYLFsOgTU4oNzcdylWX9TM7ntscPLKL9+Y9xiJjR8Fy+zo5cGvcIG7p9rDJ6cSbdHztV/adyGLSHa3p1kRbLZhpxZEVvLXmTTad2AxAuMNBv/QMbk0/RaWoRq7Rl6Y3Qli10g/jyIXtc2HFBDiwzHWbLRAuf8hVgAJKPnVge8p23t/wPgv2LwBcI0DX1buOB1o+QLXQMvjePJyKTwl5WvFJTTvGq7PuYoFlD/Yzc3ha2YMYdOFIOra53uR0pefoiUTGz3mY+cbWghGgFtkB3NV6FJ3a3GByOvF0+09k0uG1RfhZLax9pithQRpRNMP+9P28vHIsSw//DkCw00n/tAwGnLIT1vxmaH0nxLRy38hOce1bCr+Ohf1LXZ+H14Je46BBt/N62p2pO5mwfkJBAQqyBdG/SX8GNR1EqH/prLT1Bio+JeQpxcfpcPDh92P44tjXHPV3vfFfYLfSr/7/cX2nwSanKzt7D27hzblDWeKfRJ7Fgs0wuDKvBiN6TyY6SqtwpGQ+Xr6Pp7/bzCXxEcy8r63ZcXxOVm4Wk9a/z0dbppNrOPAzDPqmn+Ieu5XINnfDRfdAhSpmx3QxDNchsPlPQVqi67bWA6HbixB4fqe0WH9sPf9L+B9rk9cCUCW4CiMuHkH3uO5aBn8OKj4l5AnFZ/mGebyzYiQbA10bclXJc9I34lruufYFj53Dc75Wb/qZd38fyZqg0wBEnfmZ3OvDPxMpubunJ/Dz1qM81v0ChnTSJNOyYhgGP+3/ideWv0hyzkkALs86zROZDmpf+pBrMvF5HEoqVTlZ8MvzsOI91+cRdaDvp1Ct8Xk9bf7cpjfWvMGBjAMAXBZzGaMuGUVsxf9eqexLVHxKqDwXn2x7Fi99PoDvja3kWiz4GwZdnfGMuOEDne7hjI9+fIlPDn/GkTPL4JtnBzC8/Zu0btTO5GTiKXLynLR6bj6ZOQ6+f/AKmtYINzuSTziaeZQXlo5i0Zml6TVy83g8PZuOre7F0nYIBJWvv8f/aO8S+Ob/IP0g+IdC73ehyflPO7A77EzZOIUPN35IjjOHAGsAdze7m7ua3UWArfhL672Rik8Jldfi88uqr3h73XPsDnT9KppnB/DgZWO5tNn5HUv2RiczjvPqV3czz7KLXIuFEKeT6wMv47GbJ2Lz0+iP/Lvlu09w6wcriAwNYPWoLliLelJeKRHDMPhm+5e8vvplMpy5+BkG95xMZ1Cd6wjqPBpCo8yOWHxZKfDVnbBnkevzjk9Ch8fdMhdpf/p+XlzxIsuPLAegfuX6vHTFSzSMaHjez+3pVHxKqLwVn8ysDF6c2Z+5lp3kWSyEOZzcXKEzQ298Q4dw/sOqTfN5bdkItgW6tq1vmu3P4x3foeUFl5ucTMqzV+ZtY8Ki3VzfqgZv9G1pdhyvduTUEZ759WFWpLhWazW123nOrxb1e4yD6h5+ihBHHix8Fpa97fq8zV3Q8zVww1nbDcNg3r55jF05llR7Kn4WP+5tcS93N7u7SOcZ81YqPiVUnorPio3zeXnFo+wOcP34L8wOYmS3STqBZzHk5eXyyueD+CZ3LXara/Tn5uCOPNznbRVHOaeeb/3GliPpOk1FKZu38zueWz6GDCOPQKeTBzJzub3dc/g172veKq3SsOoD+PExwIDGveHGD107RrvBidMneGHFC/x84GcAGkU04qUrXvLZzQ9VfEqoPBQfp8PBm189xOdZv3LaaiXM4eS28B7c3/sVvVmX0MqNP/Hq8hEF5wG7OLsCz14/wyN2sJayc+jkaS5/+ResFlg9qguRFbQzuLudyjnF2F+HMzvJdaimebadlyIuJq7nWxAaaXK6UrL5G/j6XnDkQJMbXOXHDSM/4Br9mbt3Li+ufJH0nHQCbYGMuGgEfRr08bmVX+58/7b+913EXQ4m7+OeDy9navZiTlutNLLbmNhuMg/c+LpKz3m4pFl3Ph+4ihtohM0wWBV0ijtnX823iyebHU3KkYVbXTvyto6rrNJTCtYnraHPF12YnbQcq2FwX5aTae3fIK7PJ95besA1ubnvJ2D1h81fw3dDXLtBu4HFYqFnnZ58e923XF7jcuwOO8+veJ5HFj9S6IStUjwqPmXku8UfMmBOL1YFZWIzDHobDflk4Eqa17/U7GhewT8gkGcHfMGL9YZTLddJkr+FZ/e+wdNT+5CTm2N2PCkHFmxxFZ8ujbRLrjs5DSeTV73GgHkDOejIpHpeHlNCm/HAncvwb9TL7Hhlo0F3uGkKWGywfgb8NNKtT18lpArvdX6PR9s8ip/FjwX7F9Bndh/WJa9z6+v4ChWfUuZ0OBj76UDG7H2TZD8r1XMNnq8znOcHfqmTcJaCXlcM4pPec7nEXoE8i4VvrdsYOOVSdh/cbnY0MVF6di4r9pwAoGtjFR93Sc9J56Fvb+LNrR/hsECPrBy+aj2K1n1mQJCPbRXQ+Fq4/n3X9ZUTYOX7bn16q8XKgCYD+Ljnx9SsUJPDmYcZOG8gUzZNwYdmrLiFik8pOnEyifsmt+OzvDXkWSxcZA/ho95zuab9ILOjebXoqFg+uPt37gruQJDTYGNQLvfNu4Hvl35qdjQxyZIdx8h1GNSpEkqdKue34664bEtezy0zO7MofSf+hsEzzsq8cssCKra83exo5mneBzqPdl2f9wRsn+f2l2ga1ZQvr/mSHvE9cBgO3ljzBo8sfoTM3Ey3v5a3UvEpJcs3zKP/l11ZEZiBzTC42dqCD+9aRnRULbOj+QSL1cqwm8fzv+bPUz3X4Ki/ldE7x/LajPv1f0c+KP8wV1cd5nKLb9dN4vYfbyfRmU1MXh4fx/Siz4BfsVTWbsNc8bDrrO6GE2bdBcd2uP0lKgRU4JV2r/D0pU/jZ3Ud+rrth9sKzmwv/07FpxRM+u5phq15hAMBEJHn5OnY+3n6jk80gdkE7Vtfz/TrvqeVPZAcq4WPcn7jgQ86kZZx0uxoUkZyHU5+3ZYM6DDX+cp15PLsj3fx9Pp3sFvgCruDmVe8TpNur7htJZPHs1hcJzOt3Q5yTsEXd4D9VCm8jIWbL7iZaVdNo2pwVfak7eG2H25j4f6Fbn8tb6Pi40ZOh4NRU6/nnZPfknVm1dYHXT7lxivvNzuaT6tepTZT71pObxpgMQyWBJ7gzhnt2bBrjdnRpAys3ptCenYekaEBtIqtbHYcj5WSdYy7v+jKV8dWYTEMhjjCeLfvfCpd4CMTmIvD5u+a7FwhGo5tg++HuU54WgpaVGnBzGtm0rpaazJzMxm2aBhv/fEWDqejVF7PG6j4uElq2jHu+fByZlt3AdA5N5pp/X+nQVxLc4MJADabP88PmMXj0f0IczjZGWjw0OL+fL/0M7OjSSn7YeMRADo3qopNp6gokR1H1nDrl934I+cEFZxO3o24lMEDFmMN1yaQ/6hCVegz1bXSa+OXsGZaqb1UVHAUH3T7gNsbueZXfbjxQ4YtGkZWblapvaYnU/Fxg827E7hzZmdWBWXiZxgMCLycN+9eQEhQqNnR5G/6XfUk77YdT60cg+N+Vp7d+RLjZ7l36amUH3kOJ/M2JQFwdfMYk9N4pl83f8YdPw3kMHnUynXwadOhtLvWfTsUe7W4y6DzM67rP42ClNKbg+Nv9efxix9nbLuxBFgDWJS4iP5z+3Pk1JFSe01PpeJznn5c+hFDFg1gd6BBuMPJyJqDefSWiWbHkn/RqlEnpt3wIy3sAWRbLbx/ag5PTrmRvDwNDXubFXtSOJGZQ+UQfy6r68Wb6JUCwzD4cNETPLT6JbIscEmuwYyuk6hz0X1mR/Mslw2FuCsgNxO+GQylfAjq6jpXM7n7ZCKCItieup1bf7iVDcc2lOprehoVn/Pw3tcjeGbnq5zwsxKbA29f+h59ujxgdiwpgqqRsUwZuJSuedEAfG/bweDJ7TiRfsLkZOJOP2w8DMBVTavjZ9Ofu6LKybMz8tubeGv/DxgWC32dIUy4eT7hcVeYHc3zWK1w/QQICIPElfD7W6X+ki2rtmRGrxnUr1yfE9knGPTTIObunVvqr+sp9JegBJwOB89Mv5kJGXOxWy20zA5k8vU/cmHjDmZHk2IICAhm3F0LuDPoUmyGwcqgDO79/Eq27ttkdjRxg1yHk7kFh7mqm5zGc6RlneC+zzvzffoObIbBU0H1eOqOJfhX1KHCEqsUCz1ecV1fNBaO7yz1l4ypEMPHPT6mQ80O2B12RiwZwXvr3tN2Hqj4FFu2PYuhU7rwDVsB6JJbncl3/q79eTzY8L4f8Ezs3YQ5nOwIdDLs51tYsna+2bHkPC3bfYKTWblEhgZwSXyE2XE8wuETO+j/ZVcSHGmEOp28V/Ma+t78Nfhpl/nz1vI2qNfVdTLTH4aX2iqvvwr1D+WtTm8xoPEAACasn8BTvz9FrjO31F+7PFPxKYYTJ5O4Z1p7Fgccx2IY9LW24o275+vUE17ghiuH8dbF/yM618lhfwuj1j7Ml79MMTuWnIcfNrgOc/VoFq3DXEWwZf8i+s25iT3kUjXPwfSWj3JZl7GufWnk/Fks0PM18AuCvUtgw8wyeVmb1cajFz3KM22fwWaxMXv3bIb8PMSnd3rWX4Mi2rF/A3d92Y11QXYCnQYPVL6Op+74yOxY4kYXNb2KSVfNoG6OhZM2K6/uH8e7s0abHUtKIDvXwdyNrsNcvZrpEM1/WbLxIwb++gDHLQb185x82mk8F7S60+xY3iciHjqMcF3/aRRkpZTZS/dp0Ie3r3ybYL9glh9ZzsB5AzmWdazMXr88UfEpgqXrfuT+BbeyO8C1cuvpOsO497oXzY4lpSC+ZnOm9Z1PyzMrviZlzOLZ6YN0XNzDzN9ylAx7HjUqBesw13/4cunzDF3zKqctFi7NszL92q+IrnOl2bG8V9sHoUpDyDrumu9ThtrXbM+U7lOICIpgW8o2bv/xdvac9L3TXKj4/Idvfp3IiD8e46i/68zqr1/0Btd1uNvsWFKKKlWMZvLAJXTMjcBpsfAVqxn+wdXYc/LMjiZFNGvNQQBuuLAGVm1aeE6GYfD2j/fy3O4vcFgsXGuE8l7fnwmr0sjsaN7NL+DPic4JU8pkovNfNY1qyic9PiGuYhyHMw9zx9w7+OPoH2WawWwqPv9iypzneHHfeDJsVurbLbzX7XMubdbN7FhSBgICQnl70K/caG0AwM+BB7h/SntSMtJNTib/5Wh6Nr/tdA3h33ChdhY+F4cjj2dnXc8Hx5YD8H8BNXmh32L8K1QxOZmPqNMRGlwFzjxY8EyZv3ytirX4qMdHNK/SnPScdO6Zfw/z9/nOgg4Vn3/w5hcP8vaJL/5crt53IfVim5odS8qQxWplzB2z+L+wjtgMg1XBGTzwaScOJB82O5r8i2/XHsJpQJu4ysRHaff0v8vJPc1jX3RnVuZurIbBM5XbcP8tP2Lx1yKNMtX1edfpLLb/6JrsXMYigiL4sNuHdKrViRxnDo8ufpSZ28pmwrXZVHzO4bmP+zH59CIcFguX2cN5f8BiKofr/4R81f03vMNTNe4g0GmwMTiHh7/pwea9W82OJedgGAaz/nAd5rqxtUZ7/i7zdAr3f96ZBTnJ+BsGr8deQ59rp2rllhmqNIA2g1zX5z9VJsvb/y7YL5g3Or7BzQ1uxsDghZUvMHH9RK+f06ji8xdOh4MRk3vxpdO1vXfX3BjeHfSrzrkl3NT1cV5pNNy110+QkxE/38zvG5ebHUv+ZuOhNHYcPUWgn5Ve2rSwkNT0RO7+ohsrnRkEO528e8GddL2ybCfXyt90fAICKsCR9bDte1Mi2Kw2nrr0KQa3GAzAu+ve5eVVL+M0nKbkKQvFLj7bt29nzJgxdO7cmbp161K9enWaN2/OgAED+Oyzz7Db7aWRs9Tl5Nh5cMqVzPU7AMD1NOL1O3/Ez08n4hOXzpcO4q2LxhKZ5+RAAIxeeTffLzPnj5Wc26crXP9+r2oaTcUg/dvNl3RsKwNmXcMm7FRyOJl84WO0bfuI2bEkNAoucRUOfh0LTnPKhsViYUjLITxx8RMAfLbtM5787Umv3eiwyMVn7dq1dO3alRYtWrBkyRIuuugihg0bxvPPP8/tt9+OYRiMGjWKmJgYXnnlFY8qQGmnUrhvanuWBKRgNQz6B1zGcwO+wGqzmR1NypmLml3LB1d+QEyuwVF/K69ufZyP508zO5YAaadz+W79IQDuuDTO5DTlx57EZdzx/c3stTqo5nAy/YqXadZioNmxJF/bIRBYEZI3w9bZpkbp16gfY9uNxc/ix497f2ToL0M5nXfa1EylwWIU8WBeXFwcjz32GLfddhsREf+8L8by5ct54403aNmyJSNHjnRbUHdIT08nPDyctLQ0KlasCMDRE4d4aNbVbA7Mw98wuK/StdzX+yWTk0p5dzhpC/fPuYXdAQYhTid3hA/ggRtGmB3Lp01Zupfnvt9Cw+gw5j7UDovmrbB55w/839LHSbVaqO0wmNTlfarHXm52LPm7X1+Cxa+49vf5v2VgNfd/upccXMIjix4h25FNyyotGd95POGB4aZmOtf7d0kVufjk5OQQEBBQ5Ccu7v3Lwt9/cHsSNzN83q0Fb14P17ibW7o9bHZM8RBpaQcZ/MU1bArII8BpcIN/L0b2e1lvuCYwDIPO4xaz51gmz/duqhEfYOWGjxi65lWyrBaaOCy8d/UMIqo2MTuWnMvpk/Bmc7CnwU1ToekNZidibfJahiwcQkZOBvUq1eP9ru9TNaSqaXncWXyKfKirqCUmKyurWPc3y4YdyxjyU9+C3ZjHNHhCpUeKJTy8JlPv+JVLckPIsVr4Iu8HRk25hzyH904KLK+W7z7BnmOZhAbYuL5VDbPjmO7nleP4vz9cpecSpz+Tb/xRpac8C64El/6f6/rvb5mywuvvWlVtxbSrplEluAq7Tu6i/9z+7E/fb3YstyjRqq6OHTty8ODBs25fuXIlLVu2PN9M/+q9994jPj6eoKAgWrduzW+//Vbs51i+fh7DfruHg/4WquQ5efnC1+hx+R2lkFa8XVBQJSYOWEwXZyROi4U5fit57IM+nM5xmB3Np3y8wvUH+foLa1Ah0M/kNOaatWgkj2ydQq7FQhdcuzGHhmtpf7l38T3gFwxH1pmyr8+5NKjcgI96fERsWCyHTh2i/9z+bD3h+Vt5lKj4VKxYkebNm/P5558D4HQ6GTNmDO3bt+faa691a8C/mjlzJsOGDWPUqFGsXbuWdu3a0aNHDw4cOFCs53lmwzMc87NSK8fgzfaTuaJlz1JKLL7Azz+IcQN+4QZrbQB+Dt7B8A+vIvVUtrnBfMT+E5n8tNl1QtI7Lq1tbhiTTZk7mDH75+C0WLjRFsXrt/5KQIjOVeYRQqOgVT/X9WVvm5vlL2qG1WR6j+k0jGhISnYKg34aREJSgtmxzkuR5/j83cSJE3n00Ue59tpr2bdvHwcOHGDatGl06dLF3RkLXHLJJVx44YVMmDCh4LZGjRrRu3dvxo797/0o8o8RNprQiIZWP964+iviYhqUWl7xPW9+2Z/JWWsBuCgznGdvnkOtqMomp/JuT3+7iY9X7KfjBVWYdufFZscxheF0Mu67W5iW7vq/8buC43noxq+x2Hx79MvjpOyBd1qD4YTBv0N0+TlbQEZOBg/+8iBrjq4h0BbIa+1fo1NspzJ7fVPm+Pzd4MGDefDBB/n8889JSEjgiy++KNXSk5OTw5o1a+jWrfC5srp168ayZcvO+Ri73U56enqhD4Am2f6832eeSo+43bA+H/FIZBeshsHq0DSe+KIrm/YfMjuW10rJzOHLNYkA3Nu+jslpzJGXl8MzX/YsKD2PVr6QYX2+U+nxRBF1oPF1ruvL3jE3y9+EBYQxsctEOtbqiN1h5+FFD/Pdru/MjlUiJSo+qamp3HjjjUyYMIH333+fm2++mW7duvHee++5O1+B48eP43A4qFatWqHbq1WrRlJS0jkfM3bsWMLDwws+atWqBcCbfedQpXJMqWUV3zbw6jd4LvZWApwGG0LtPDu3F79t2mJ2LK/00fJ9ZOc6aVYjnLZ1Is2OU+bs9gwe+bwz32YfwmYYPF+9CwOuna5TUHiyy4a6LjfNglPJ5mb5myC/IN7o+AbX1r0Wh+Hgqd+fYvrm6WbHKrYSFZ+mTZty9OhR1q5dyz333MMnn3zC5MmTefrpp+nVq5e7Mxby96XChmH84/LhJ598krS0tIKPxETX/xmGhVYq1Ywi1105inGNHyTE6WRbsIPXlvXl2+W/mx3Lq5zOcfDRctek5nvb1/G5bQROZSbzf5935hfHSQIMg3F1b6F3tzfMjiXnq8aFUKMNOHPhj/JXKvysfjx/+fP0b9wfgNcTXuetP97yqPN7laj4DB48mCVLlhAfH19wW9++fVm/fj05OTluC/dXUVFR2Gy2s0Z3kpOTzxoFyhcYGEjFihULfYiUlQ6X3Meki5+nksPJ3kCYsOk+Js+fY3Ysr/Hpyv2kZOZQs3IwPZpGmx2nTJ1I2c2gL7qzmtOEOp1MbPYAV7Z7yuxY4i4X3+O6TJgGjjxTo5yL1WLl0TaP8tCFDwHw4cYPeW7FczicnrGatUTF5+mnn8ZqPfuhNWvWZMGCBecd6lwCAgJo3br1Wc+/YMECLrvsslJ5TZHz1aLJDUzv9B7V8gwOB1j4JPEJXvtqqkf931F5lJWTx8TFuwF4oFM9/Gy+c77lQ4fXMOC769lqzSPC4WTKJc9yUevBZscSd2rcG0IiIf0g7Jhndppzslgs3N3sbka3HY3VYuWrHV/x2JLHyHGUzuCHOxX5r0Vxl4wfOuT+CZ3Dhw/nww8/ZMqUKWzdupWHH36YAwcOMHiw/tFL+VUnvgMf9/qY2nlw3M/Kt+mv8/RH/yNXGx2W2EfL93P8VA6xESHc2Np39qjZtXch/X8awH6rQYzD4KNO42nc+CazY4m7+QdBqzN7y63+0Nws/+GmBjfxeofX8bf6s2D/Au5feD+ZuZlmx/pXRS4+F110Effccw+rVq36x/ukpaXxwQcf0LRpU77++mu3BPyrvn378uabb/Lcc8/RsmVLlixZwo8//khcnLanl/KtenQrPr5hNo3ybKTbrMx3TmXEB0+RlVP+hrHLu1P2PN4/M9rzUOf6+PvIaM+6zV8wYNFDJFst1HNY+KjHR8TFl91yYiljbe4ELLDnVzi+y+w0/6prXFfe7fwuwX7BrDyykrt/upvU7FSzY/2jIu/jk5KSwksvvcSUKVPw9/enTZs2xMTEEBQURGpqKlu2bGHz5s20adOGp556ih49epR29mJz5z4AIiWRdSqZIV/0JMFmx88waJvRjufvfJPICoFmR/MY7yzcyf8W7KBOVCjzH27vE4e5lia8x/CN73HaaqGF0493e88ivLJvLt/3KZ/eDDt/gssfgq7PmZ3mP208tpH7F97PSftJ4sPjmdR1EtGh7pl/Z8o+PgcPHuSVV17h8OHDTJw4kQYNGnD8+HF27twJQL9+/VizZg2///57uSw9IuVBSIWqvN/vVzoRTp7FwtKw33jywwHsOXbK7GgeITkju2Buz0Nd6vtE6flxyXM8uMlVei4nmEk3/6TS4yta3e66XP95uZzk/HfNqjRj+lXTqRZSjb1pe7lj7h3sSdtjdqyzFHnEJ39FVZUqVahTpw6rV68mMtKz9s3QiI+UF468HMbMvJpv844AcFFqLPdcPY229aqYnKx8G/HVer5IOEiLWpX45v8uw2r17iXsn/00lJeP/IJhsdDDWokXb56Lf2AFs2NJWcnLgXENIesE3PYFNOhudqIiOXLqCPcuuJd96fuoHFiZCV0m0CTq/E6Sa8qIT6VKldizx9Xc9u3bh9OpiZkiJWXzC+C5W+cxKLQRAKsrH2DiD334cvU+c4OVY5sOpfHlGtfJkZ+5urFXlx7D6eS97/oxNulXDIuFWwNr8vJtv6j0+Bq/AGh2s+v6uk/NzVIM1StUZ3qP6TSJbEKqPZVBPw1i5ZGVZscqUOTic+ONN9KhQwfi4+OxWCy0adOGOnXqnPNDRP6bxWrl4Zu+4NEo13YMf1Q6wbcr+vL6j+txOrXc/a8Mw+C5OVswDLiuZQyt47z3/GeOvBxe/OoaJpzcAMD9FZvw5M0/YLX5m5xMTJF/4tLtcyErxdwsxRARFMHk7pO5JPoSsvKy+L+f/4+f9/9sdiygmCcpnTdvHrt27WLo0KE899xzhIWFnfN+Dz30kNsCupMOdUl59d0vIxl9YDYOi4XGmf5UCxvHq7e0I8jfZna0cmHWmoM88uV6gvyt/PJIR2IqBZsdqVScPp3KiK+uZpEzHYth8GS1dtzaY8J/P1C828R2kLQBerwKl9xndppisTvsPLHkCX4+8DNWi5VnLn2GGxvcWOzncef7d4nOzn7nnXfy9ttv/2PxKa9UfKQ8W7xiHI9snYLdaqFutoWQvKd4q/91VAnz7RVfx0/Z6TJuMSezcnms+wUM6VTP7EilIjV1Dw98dxMbLLkEGAYv17mZru2fMTuWlAcr34e5I6B6C7hvidlpis3hdPD8iueZtXMWAINbDOb+FvcX6zQzpp+dferUqR5XekTKuw6XDmfShSOo6HSyO8jgeOBz3P3eRNYlnjQ7mqnGzN7MyaxcGlev6LVnYE88uJw7vunNBksuFZ1OPmz5qEqP/KnpTWCxwZH15X5Pn3OxWW2Mbjuae5q5TsUxcf1Envr9KXIduabk8f61oCIe5MIW/fm44zvEOAyO+Fs4UWUST01/lS9WJ5odzRQLthzl+w1HsFktvHpTc6/crHDT1q+5fcE97Le5dmP+uMNbtGo50OxYUp6ERkLdM5tVbnb/5sBlwWKxMPTCoYxuOxqbxcbs3bO5f+H9ZORklHkW7/srIuLh6sRfyafXfkVjp400m5WkGrP5YsGzPP3tJnLyfGc1ZXJ6Nk/Mck3wvaddHZrWCDc5kfstWfU2g1Y8Q4rVQiOnlU+unkGdOl3MjiXlUdMz82I2fgUefK6/mxrcxDtXvkOwXzArjqxgwLwBJGUm/fcD3UjFR6QciopqyNS+P9PeUgG71cKOmBXs3fkU/T5cwbEMu9nxSp3DaTBs5jpOZObQqHpFhnWpb3Ykt/vy50cYumUSp60WLiOYqTfNo0rVZmbHkvKqYS+wBcDx7ZC8xew056VdzXZMu2oaUcFR7EzdSb8f+7E9ZXuZvb6Kj0g5FRISxVu3/krfwJoYFgvrq+4mOPsJrn7rF5btPm52vFI1cfFulu0+QUiAjfG3tfKq1W2OvBxe/fI6njs0H4fFwrV+VRh/66+EhlU3O5qUZ0HhUL+b6/qmWeZmcYPGkY35tOen1A2vS3JWMgPmDWDZoWVl8toqPiLlmJ9/EKNu/oHhkRcDsLZyKrUjRnHflHmMW7ADhxfu9/PbzmOMW7ADgOeua0rdKt6zad+pU0kM/awDH2e5NoN9ILw5L9z6M/4BoSYnE4/Q9AbX5aavPfpwV76YCjF81PMjLoq+iMzcTO5feD8zts2gBIvNi0XFR6Scs1it3Hn1ZF6L70Og02BLaC41a7/MrN++59YPVpCUlm12RLfZlZzB/Z/+gcNpcFPrmtx4YQ2zI7nNocMJ3PFld5YYpwh0GrwefxP39f4Ui1V/hqWIGlwF/iGQuhcO/2F2GreoGFCRiV0mcm3da3EYDl5a+RLPr3ieXGfprfjSvzgRD3FV+2eYdsloqjoMDgZYoPZkTid/SY+3ljBv0xGz4523lMwcBk1LICM7j4tqV+bF65sWa5+P8mzdphnc9tNAdlmdRDkMpl3yDN3bjzY7lniagNA/z9e19Xtzs7hRgC2AFy5/gUdaP4IFC1/u+JL7FtzHyeyTpfJ6Kj4iHqRp4z7M6PU5TZ1+pNusHKw1n8ZBkxj8yR8M+3wtaVnm7ItxvjKyc7lz2moOpGRRKyKYibe3JtDPO+b1zP51FHclvEiK1UJDp5UZPT6maeObzY4lnqrh1a7LbT+Ym8PNLBYLA5sO5J0r3yHUP5TVSau59Ydb2ZXq/n2LVHxEPEzVak2Zessv9PSLJM9iYUO17bSv/jLfr9tL1zcW8+u2ZLMjFktWTh53TUtgfeJJKoX4M2XARURW8PzdqnNzsnhxZi9GHZhNjsXClZYwpvdZQHT1VmZHE09WvytY/V2ru47vNDuN23Wo1YFPenxCzQo1OXjqILfPvZ3FiYvd+hoqPiIeKCi4Mi/f+gsPRbTGYhisrXSSFvHPYpzewZ3TVvPgjLUcTS//c3/STucycMpqVu1LISzIj0/uuoT61Tx/V/jk5M0M+qwdn2cfAOD/KjbhjX6LCalQ1eRk4vGCwiG+veu6l4365KtXuR6f9fqMNtXakJmbyYO/PMiHGz902/Or+Ih4KIvVyt3XTOOtCwZSwWmwI8hJcJ13aR66kDnrD9P5f4uZsnQveY7yuelhcno2fd9f7io9gX5Mu/Nir9ik8I+Nn9L3+76ss+QQ5jQYf8FA7r/+c51dXdynYS/XpZcWH4DKQZWZ1HUSfS/oi4Hh1uJTopOUeiqdpFS81YEDS3l44RB2WJ1YDYMrMuvyY+IgDPxoGB3GqF6NaFe/itkxC6w9kMrgT9ZwNN1OVIVAPhp0MY1jPPvfpOF0MmP+UF5LWkSexUI9h4U3O48nLq692dHE26QfgXENXdcf2Q5h0ebmKWXf7PyG91e9z0+3/2Te2dk9lYqPeLPTWSm88G0fZue65vhcmhfClkMPcSjLNYpyRb0oHr+qIc1qmjeq4nQafLR8Hy/9uI0ch5N6VSsweUAb4iI9ex+btLQDjJ59GwudaQB0t1Xmud5f6dCWlJ4ProRDa+DqN6DNILPTlLoTJ08QVTnKvLOzi0j5ExwSwQu3LGB0TFcCDIMVflmE1nyJu5ptxt9mYemu41wzfin3fZzAHwdSyzxfYkoW/aesYsycLeQ4nHRvUo1vh1zu8aVn7cZP6TOrFwudafgZBiOqXMZrty1S6ZHS5QOHu/7K3+q+Q8Ua8RHxQpu3fc0jy0ZzyAZWw2BAaGMS8x7h2w3HCzZ8vbh2BPe0r0OnC6rgV4pnPc+05/H+kj28v3g39jwnQf5WRvZsxO2XxGG1eu4+PY68XKb8eDfvpqzBYbFQywGvtR1Dk0Y3mh1NfEHyVnjvUrAFwuP7ICDE7ESlyp3v3yo+Il7qVMYRxs7px+zcYwA0d/rzYJvX+Hp7FN+sPUSuw/VPv1rFQG68sCY3tq7p1tNDnDhlZ/qyfUxfvp+00679hdrWieSF6z3/NBSHDyfw9IL/YxWulXM9/SJ5+toZVND5tqSsGAa80RTSD8JtX0KDbmYnKlUqPiWk4iO+aO7iMTy/5ysyrBZCnQaPx/bk0gvHMG35fr5ISCT1L5se1q0SSpdG1eh4QVVa1AonJMCvWK+VmpnDb7uOM3vdIRbvOFZQruKjQhnR/QKuahrt0bsxG04nX/8ygtcS55FptRDsNBgZ24vrOo3VqSek7M0eCn9Mh4vvg56vmp2mVKn4lJCKj/iqw4cTeHL+ffxhyQHgckJ4puu7RFZtxS9bk/kiIZHfdh4n7y8nPbVZLTSoFkbTmIrERoRQMyKY8GB/QgL8sABZuQ4ysvNITMli3/FMNhxMY/vRjEKv26JmOPd1qEv3JtHYPPiwFsDRoxsZ89M9LDUyAWhp+PNCp7eIi2tncjLxWVvnwMzbIaIODF1rdppSpeJTQio+4ssceTlMnzuY946vwm61EOI0GFa9I327vonV5kfa6VyW7DjGwq1HWbEnhaQSboBYr2oFrmoSzXUtY7xiM0KnI49vfnmc/x38iQyrhQDD4MGoS7njqvew+QWYHU98WXY6vBoPzjx48A+IrGt2olKj4lNCKj4isG/fYkYvGl4w+nOhEcCoK16gQb0ehe6XlJbNusRUdhw9RWJKFodOnuaUPY9Mex6GASGBNkIC/KhZKZi4yFAuiK5Am9oRRHnB6Sby7dg1l+eXjmKdxXU4sInTxovtX6Nu3a4mJxM5Y9rVsO836PEqXHKf2WlKjYpPCan4iLg4HXnMXDCMN44s4rTVgs0w6BtSmyHdJ1AxvJbZ8UyXlXmMCXPv5eNTO3FYXHN5hlS9jH7d3sbPP8jseCJ/Wvom/Dwa6neDfl+anabUuPP9W7PxRHyQ1ebHrVeN57vu0+lqDcdhsfDZ6f1c83UPvv75UZyOPLMjmsKRl8M3Cx/nmpmdmJa5C4fFQhdrOLOv+pgBvSap9Ej5U6+L63Lvb5Bb/s/PVx5oxEdEWL5mIi+vf489NtefgwucVoY2HkS7ix70mdVKyxLe438b32eH1XVusxoOGNlkEO0vedjkZCL/wjBgXGPIOAy3fw31OpudqFRoxEdE3Kpt68F8dccKHo1qSwWnwXarkyHbPmTgRxfxx/qPzI5Xqv5Y/xH3TL+Y+zZPYIfVSZjT4NGoS5ndb5lKj5R/FgvU6ei6vneJqVE8hUZ8RKSQk6l7mbxwODMydmI/swT9EoIY1ORO2l442GtGgFavnczE9RNYZbED4GcY3BJah/u6vEWlyvEmpxMphvWfwzf3QcyFcO+vZqcpFZrcXEIqPiJFd/ToBt7/dQTfZB8k78ymg42cNgbVvZ4ubR/3yPkuublZ/LzsVT7dM5v1VtdKLT/DoHdQDe664llq1rzU5IQiJZB+GMY1AosVRuyF4EpmJ3I7FZ8SUvERKb7DhxP4aOmzfJ25l9NnRoCqOgyur9yUGy55lJiYNiYn/G/Hj23lq2Uv8uWJdSTbXN+Dv2FwQ1At7mr/AtVjWpucUOQ8vdMGTuyEWz778wSmXkTFp4RUfERKLjVlN58veZoZqRtIPVOALIbBZZZQronrTsc29xNaIdrklH86nZXCL6veZM7+n1huZOI8M2oV5TC4OaI5fS5/iqgqjU1OKeIm3w+HhMlwyWDo8YrZadzOJ4vPiy++yA8//MC6desICAjg5MmTxX4OFR+R85djz+CXlW/w1d45rOTP5bMBhsEVtnC61OxI22Z3EBXVsMyzpaTs4rd1H7Lo0FJ+zztZMEIF0MLpz63xV9Ot7eP4B4aWeTaRUrX5W/hyAFRpBENWmJ3G7Xyy+IwePZpKlSpx8OBBJk+erOIjUg4cOLCU79a8w/zULeyzFf7aBU4rbcPq0LL6JTSt15Nq1Zq7/fVTUnaxdtss1hxaxppT+9lmySsY2QHXkvRrKjfl6gvv1zm1xLtlpbhOXwHw6E6oUNXcPG7mk8Un37Rp0xg2bJiKj0g5Yjid7Ng9l582Tmdp2g62Wh1n3aeKw+ACWwVqB1chrmJtYiMbEhkeR+Xw2lSuFH/WKIzhdGK3p5GZeYzjqTs5mrKDpJP7SMw4wM7Mw+xwZHLMdvaJTxs5bXSs1IiOjfrSqMG1XrMKTeQ/TbwCkjbCjZOh2U1mp3Erd75/+7kpU7lkt9ux2+0Fn6enp5uYRsR7WaxWLqjfiwvq92IorpGYFRs/ZtXh5WzMPsoui4NjNgvHyGTp6Uw4vQ+OLir0HP6GgZ8BNnCd/d0CDsu/nNH9TOmp67DQOiSG1tUvoXXDG0tlZEnEI8R3cBWfvYu9rvi4k1cXn7Fjx/Lss8+aHUPE50RE1KNnh2fpeebzrKzjbNv5I7uPrmV/2l72n07moCOTFByknSk4uRYLuf/Qcyo5DaKxEW0LJSYognqV6tMg5hLq1e5YriZUi5gqvj0sHw/7fjc7SblmavEZM2bMfxaT1atX06ZNyZbLPvnkkwwfPrzg8/T0dGrV0gkYRcpaSEgUF7boz4X0P+trTkceGRkHOZ19EocjB4cjB6fTQUhwJCEhkQQHR2LzCzAhtYiHqXUxYIGU3XAq2evm+biLqcXngQce4JZbbvnX+9SuXbvEzx8YGEhgYGCJHy8ipc9q8yO8Um3CzQ4i4umCK0PVxpC8GQ4sh8bXmZ2oXDK1+ERFRREVFWVmBBEREe8Re+mZ4rNCxecfeMxyhwMHDrBu3ToOHDiAw+Fg3bp1rFu3jlOnTpkdTUREpHyIu8x1eWC5uTnKMY+Z3PzMM88wffr0gs9btWoFwK+//krHjh1NSiUiIlKOxJ4539yRDWA/BYEVzM1TDnnMiM+0adMwDOOsD5UeERGRM8JrQngtMBxwcLXZacoljyk+IiIiUgT5oz4HvO/UFe6g4iMiIuJNYtu6LjXP55xUfERERLxJfvE5uBocueZmKYdUfERERLxJlYYQVAlysyBpg9lpyh0VHxEREW9itULNM2c8OLjG3CzlkIqPiIiIt6lxpvgcSjA3Rzmk4iMiIuJtCkZ8VHz+TsVHRETE29Ro7bpM2Q1ZKeZmKWdUfERERLxNSARE1HFdP/SHuVnKGRUfERERb6R5Puek4iMiIuKNNM/nnFR8REREvFHBiM8aMAxzs5QjKj4iIiLeKLop2ALgdAqk7DE7Tbmh4iMiIuKN/AIhurnruiY4F1DxERER8Vb5y9o1wbmAio+IiIi3qvmXeT4CqPiIiIh4r+otXZdJm8DpMDVKeaHiIyIi4q0i64J/KOSdhuM7zU5TLqj4iIiIeCurDaKbua4fWW9ulnJCxUdERMSbVW/hulTxAVR8REREvJuKTyEqPiIiIt4sv/gkbQCn09ws5YCKj4iIiDercgHYAsGeDql7zU5jOhUfERERb2bzh2pNXNd1uEvFR0RExOtpnk8BFR8RERFvp+JTQMVHRETE2/21+BiGuVlMpuIjIiLi7ao2BqsfnE6BtESz05hKxUdERMTb+QdB1AWu60c3m5vFZCo+IiIiviB/ZZeKj4iIiHi9ao1dlyo+IiIi4vWqNXVdJm8xN4fJVHxERER8QdUzIz7Hd0Ke3dwsJlLxERER8QUVYyCoEhgOOLbd7DSmUfERERHxBRaLJjjjIcVn37593HXXXcTHxxMcHEzdunUZPXo0OTk5ZkcTERHxHPnFJ9l3i4+f2QGKYtu2bTidTt5//33q1avHpk2buOeee8jMzOT11183O56IiIhnqKqVXR5RfK666iquuuqqgs/r1KnD9u3bmTBhgoqPiIhIUeWv7Drquyu7PKL4nEtaWhoRERH/eh+73Y7d/ufM9fT09NKOJSIiUn5Vbei6PJUEmScgNNLcPCbwiDk+f7d7927eeecdBg8e/K/3Gzt2LOHh4QUftWrVKqOEIiIi5VBgGFSKc1330Xk+phafMWPGYLFY/vUjISGh0GMOHz7MVVddRZ8+fbj77rv/9fmffPJJ0tLSCj4SE337xGwiIiK+frjL1ENdDzzwALfccsu/3qd27doF1w8fPkynTp1o27YtkyZN+s/nDwwMJDAw8HxjioiIeI9qjWH7D3B0k9lJTGFq8YmKiiIqKqpI9z106BCdOnWidevWTJ06FavVI4/SiYiImKtqI9elj25i6BGTmw8fPkzHjh2JjY3l9ddf59ixYwVfi46ONjGZiIiIh6lyZoLzse1gGK6NDX2IRxSf+fPns2vXLnbt2kXNmjULfc0wDJNSiYiIeKDIemCxgj0NTh2FMN8aQPCI40UDBw7EMIxzfoiIiEgx+AVC5XjX9WPbzM1iAo8oPiIiIuJGVS5wXR7bYW4OE6j4iIiI+JqC4qMRHxEREfF2UWeKz3GN+IiIiIi3Kxjx8b0l7So+IiIiviaqgesyMxmyUszNUsZUfERERHxNYAWoeGZ7GB873KXiIyIi4ot8dIKzio+IiIgv8tEl7So+IiIivkgjPiIiIuIzfHRJu4qPiIiIL8of8UlLBPspc7OUIRUfERERXxQSAaFVXNdP7DQ3SxlS8REREfFVkfVclyd2m5ujDKn4iIiI+KqIuq5LFR8RERHxepFnik+Kio+IiIh4u4JDXbvMzVGGVHxERER8Vf6Iz4ldYBjmZikjKj4iIiK+KqKO6zI7zWdOVqriIyIi4qv8gyG8luu6jxzuUvERERHxZfmjPio+IiIi4vXyJzj7yMouFR8RERFf5mMru1R8REREfFmkb21iqOIjIiLiywoOde0Bp9PcLGVAxUdERMSXVYoFqx/kZkHGEbPTlDoVHxEREV9m84dKca7rPjDPR8VHRETE1/nQyi4VHxEREV9XsLJLxUdERES8XUS86zJlj7k5yoCKj4iIiK/LLz6p+0yNURZUfERERHxd5b8UHy8/S7uKj4iIiK8LrwUWq2tJ+6lks9OUKhUfERERX+cXABVruq57+eEuFR8RERGBymf28lHxEREREa9XubbrUsWnfLj22muJjY0lKCiI6tWrc8cdd3D48GGzY4mIiHiHguKz19QYpc1jik+nTp344osv2L59O7NmzWL37t3cdNNNZscSERHxDj6ypN3P7ABF9fDDDxdcj4uL44knnqB3797k5ubi7+9vYjIREREvoENd5VdKSgqffvopl112mUqPiIiIO+Tv5ZNxBHJPm5ulFHlU8Xn88ccJDQ0lMjKSAwcO8N133/3r/e12O+np6YU+RERE5ByCK0NgRdf11P3mZilFphafMWPGYLFY/vUjISGh4P6PPfYYa9euZf78+dhsNvr374/xLztMjh07lvDw8IKPWrVqlcW3JSIi4nksFp843GUx/q05lLLjx49z/Pjxf71P7dq1CQoKOuv2gwcPUqtWLZYtW0bbtm3P+Vi73Y7dbi/4PD09nVq1apGWlkbFihXPL7yIiIi3mXkHbJ0NV70Clw42O02B9PR0wsPD3fL+berk5qioKKKiokr02Py+9tdi83eBgYEEBgaW6PlFRER8jg+M+HjEqq5Vq1axatUqrrjiCipXrsyePXt45plnqFu37j+O9oiIiEgx+cBePh4xuTk4OJivv/6azp07c8EFFzBo0CCaNm3K4sWLNaIjIiLiLj6wl49HjPg0a9aMX375xewYIiIi3u2vh7oMwzXh2ct4xIiPiIiIlIHwWmCxQV42ZCSZnaZUqPiIiIiIi80fKsa4rqclmpullKj4iIiIyJ/Cz+x5d/KAuTlKiYqPiIiI/KnSmeKjER8RERHxepViXZca8RERERGvV3CoSyM+IiIi4u10qEtERER8Rnj+oa5E114+XkbFR0RERP4UXtN1mZsJWSnmZikFKj4iIiLyJ/8gqFDNdT3N+yY4q/iIiIhIYV48wVnFR0RERArz4gnOKj4iIiJSmEZ8RERExGd48SaGKj4iIiJSWH7x0eRmERER8Xo61CUiIiI+I39yc/ZJsGeYGsXdVHxERESksMAwCKrkuu5loz4qPiIiInK2/FEfL5vgrOIjIiIiZ6sU57r0sr18VHxERETkbOEa8RERERFfkX+y0vRD5uZwMxUfEREROVt4DddlmoqPiIiIeLuKZ4pP+mFzc7iZio+IiIicrWKM6zLjMDid5mZxIxUfEREROVuFaLBYwZkHmclmp3EbFR8RERE5m83PVX7AqyY4q/iIiIjIueUf7vKieT4qPiIiInJuXriyS8VHREREzq1gZZeKj4iIiHg7HeoSERERn6ERHxEREfEZKj4iIiLiMwoOdR3xmk0MVXxERETk3MLyNzHMhcxjZqdxCxUfEREROTebP1So5rruJYe7PK742O12WrZsicViYd26dWbHERER8W5etrLL44rPiBEjiImJMTuGiIiIb/CyCc4eVXzmzp3L/Pnzef31182OIiIi4hu8rPj4mR2gqI4ePco999zDt99+S0hISJEeY7fbsdvtBZ+npaUBkJ6eXioZRUREvI4tAuwGHNkLJr1/5r9vG4Zx3s/lEcXHMAwGDhzI4MGDadOmDfv27SvS48aOHcuzzz571u21atVyc0IRERFv98mZD/OcOHGC8PDw83oOi+GO+lRCY8aMOWcx+avVq1ezbNkyZs6cyZIlS7DZbOzbt4/4+HjWrl1Ly5Yt//Gxfx/xOXnyJHFxcRw4cOC8f3ByftLT06lVqxaJiYlUrFjR7Dg+Tb+L8kO/i/JFv4/yIy0tjdjYWFJTU6lUqdJ5PZepxef48eMcP378X+9Tu3ZtbrnlFubMmYPFYim43eFwYLPZ6NevH9OnTy/S66WnpxMeHk5aWpr+IzaZfhflh34X5Yd+F+WLfh/lhzt/F6Ye6oqKiiIqKuo/7/f222/zwgsvFHx++PBhunfvzsyZM7nkkktKM6KIiIh4EY+Y4xMbG1vo8woVKgBQt25datasaUYkERER8UAetZz9fAUGBjJ69GgCAwPNjuLz9LsoP/S7KD/0uyhf9PsoP9z5uzB1jo+IiIhIWfKpER8RERHxbSo+IiIi4jNUfERERMRnqPiIiIiIz/CZ4vPee+8RHx9PUFAQrVu35rfffjM7kk8aO3YsF110EWFhYVStWpXevXuzfft2s2P5vLFjx2KxWBg2bJjZUXzWoUOHuP3224mMjCQkJISWLVuyZs0as2P5nLy8PJ566ini4+MJDg6mTp06PPfcczidTrOjeb0lS5ZwzTXXEBMTg8Vi4dtvvy30dcMwGDNmDDExMQQHB9OxY0c2b95c7NfxieIzc+ZMhg0bxqhRo1i7di3t2rWjR48eHDhwwOxoPmfx4sUMGTKEFStWsGDBAvLy8ujWrRuZmZlmR/NZq1evZtKkSTRv3tzsKD4rNTWVyy+/HH9/f+bOncuWLVv43//+d95b80vxvfLKK0ycOJHx48ezdetWXn31VV577TXeeecds6N5vczMTFq0aMH48ePP+fVXX32VcePGMX78eFavXk10dDRdu3YlIyOjeC9k+ICLL77YGDx4cKHbGjZsaDzxxBMmJZJ8ycnJBmAsXrzY7Cg+KSMjw6hfv76xYMECo0OHDsZDDz1kdiSf9PjjjxtXXHGF2THEMIxevXoZgwYNKnTbDTfcYNx+++0mJfJNgPHNN98UfO50Oo3o6Gjj5ZdfLrgtOzvbCA8PNyZOnFis5/b6EZ+cnBzWrFlDt27dCt3erVs3li1bZlIqyZeWlgZARESEyUl805AhQ+jVqxddunQxO4pPmz17Nm3atKFPnz5UrVqVVq1a8cEHH5gdyyddccUVLFy4kB07dgCwfv16li5dSs+ePU1O5tv27t1LUlJSoffywMBAOnToUOz3co84ZcX5OH78OA6Hg2rVqhW6vVq1aiQlJZmUSsB1vHb48OFcccUVNG3a1Ow4Pufzzz/njz/+YPXq1WZH8Xl79uxhwoQJDB8+nJEjR7Jq1SqGDh1KYGAg/fv3NzueT3n88cdJS0ujYcOG2Gw2HA4HL774IrfeeqvZ0Xxa/vv1ud7L9+/fX6zn8vrik++vZ3YH15vu32+TsvXAAw+wYcMGli5danYUn5OYmMhDDz3E/PnzCQoKMjuOz3M6nbRp04aXXnoJgFatWrF582YmTJig4lPGZs6cySeffMJnn31GkyZNWLduHcOGDSMmJoYBAwaYHc/nueO93OuLT1RUFDab7azRneTk5LOao5SdBx98kNmzZ7NkyRKdaNYEa9asITk5mdatWxfc5nA4WLJkCePHj8dut2Oz2UxM6FuqV69O48aNC93WqFEjZs2aZVIi3/XYY4/xxBNPcMsttwDQrFkz9u/fz9ixY1V8TBQdHQ24Rn6qV69ecHtJ3su9fo5PQEAArVu3ZsGCBYVuX7BgAZdddplJqXyXYRg88MADfP311/zyyy/Ex8ebHcknde7cmY0bN7Ju3bqCjzZt2tCvXz/WrVun0lPGLr/88rO2ddixYwdxcXEmJfJdWVlZWK2F3xptNpuWs5ssPj6e6OjoQu/lOTk5LF68uNjv5V4/4gMwfPhw7rjjDtq0aUPbtm2ZNGkSBw4cYPDgwWZH8zlDhgzhs88+47vvviMsLKxgJC48PJzg4GCT0/mOsLCws+ZVhYaGEhkZqflWJnj44Ye57LLLeOmll7j55ptZtWoVkyZNYtKkSWZH8znXXHMNL774IrGxsTRp0oS1a9cybtw4Bg0aZHY0r3fq1Cl27dpV8PnevXtZt24dERERxMbGMmzYMF566SXq169P/fr1eemllwgJCeG2224r3gu5Y9mZJ3j33XeNuLg4IyAgwLjwwgu1fNokwDk/pk6danY0n6fl7OaaM2eO0bRpUyMwMNBo2LChMWnSJLMj+aT09HTjoYceMmJjY42goCCjTp06xqhRowy73W52NK/366+/nvP9YcCAAYZhuJa0jx492oiOjjYCAwON9u3bGxs3biz261gMwzDc0dREREREyjuvn+MjIiIikk/FR0RERHyGio+IiIj4DBUfERER8RkqPiIiIuIzVHxERETEZ6j4iIiIiM9Q8RERERGfoeIjIiIiPkPFR0RERHyGio+IeKxjx44RHR3NSy+9VHDbypUrCQgIYP78+SYmE5HySufqEhGP9uOPP9K7d2+WLVtGw4YNadWqFb169eLNN980O5qIlEMqPiLi8YYMGcLPP//MRRddxPr161m9ejVBQUFmxxKRckjFR0Q83unTp2natCmJiYkkJCTQvHlzsyOJSDmlOT4i4vH27NnD4cOHcTqd7N+/3+w4IlKOacRHRDxaTk4OF198MS1btqRhw4aMGzeOjRs3Uq1aNbOjiUg5pOIjIh7tscce46uvvmL9+vVUqFCBTp06ERYWxvfff292NBEph3SoS0Q81qJFi3jzzTf5+OOPqVixIlarlY8//pilS5cyYcIEs+OJSDmkER8RERHxGRrxEREREZ+h4iMiIiI+Q8VHREREfIaKj4iIiPgMFR8RERHxGSo+IiIi4jNUfERERMRnqPiIiIiIz1DxEREREZ+h4iMiIiI+Q8VHREREfIaKj4iIiPiM/weVUOhOEx5HNwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "pylab.plot(xx, y10, label='O(10)')\n",
    "pylab.plot(xx, y20, label='O(20)')\n",
    "pylab.plot(xx, y, label='sin(x)')\n",
    "\n",
    "pylab.axis([0, 10, -4, 4])\n",
    "pylab.xlabel('x')\n",
    "pylab.ylabel('f(x)')\n",
    "\n",
    "pylab.legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Linear equations and matrix inversion\n",
    "\n",
    "SymPy has a `Matrix` class and associated functions that allow the symbolic solution of systems of linear equations (and, of course, we can obtain numerical answers with `subs()` and `evalf()`). We shall consider the example of the following simple pair of linear equations:\n",
    "\n",
    "$$\\begin{aligned}\n",
    "3x + 7y&= 12z\\\\\n",
    "4x - 2y&= 5z\\end{aligned}$$\n",
    "\n",
    " We may write this system in the form $A\\vec{x}=\\vec{b}$ (multiply *A* by $\\vec{x}$ if you want to verify that we recover the original equations), where\n",
    "\n",
    "$$A=\\left(\\begin{array}{cc}\n",
    "3 & 7\\\\\n",
    "4 & -2 \\end{array} \\right),\\qquad\n",
    "\\vec{x}=\\left(\\begin{array}{c}\n",
    "x\\\\\n",
    "y \\end{array}\\right),\\qquad\n",
    "\\vec{b}=\\left( \\begin{array}{c}\n",
    "12z\\\\\n",
    "5z \\end{array}\\right).$$\n",
    "\n",
    "Here we included a symbol, *z*, on the right-hand side to demonstrate that symbols will be propagated into the solution. In many cases we would have *z* = 1, but there may still be benefit to using SymPy over a numerical solver even when the solution contains no symbols because of its ability to return exact fractions rather than approximate `float`s.\n",
    "\n",
    "One strategy to solve for $\\vec{x}$ is to invert the matrix *A* and pre-multiply, *i.e.* $A^{-1}A\\vec{x}=\\vec{x}=A^{-1}\\vec{b}$. SymPy’s `Matrix` class has an `inv()` method that allows us to find the inverse, and `*` performs matrix multiplication for us, when appropriate:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 78,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}3 & 7\\\\4 & -2\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "⎡3  7 ⎤\n",
       "⎢     ⎥\n",
       "⎣4  -2⎦"
      ]
     },
     "execution_count": 78,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from sympy import symbols,Matrix\n",
    "x, y, z = symbols('x,y,z')\n",
    "A = Matrix(([3, 7], [4, -2]))\n",
    "A"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 79,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}\\frac{1}{17} & \\frac{7}{34}\\\\\\frac{2}{17} & - \\frac{3}{34}\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "⎡1/17  7/34 ⎤\n",
       "⎢           ⎥\n",
       "⎣2/17  -3/34⎦"
      ]
     },
     "execution_count": 79,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A.inv()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 80,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}12 z\\\\5 z\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "⎡12⋅z⎤\n",
       "⎢    ⎥\n",
       "⎣5⋅z ⎦"
      ]
     },
     "execution_count": 80,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "b = Matrix(( 12*z,5*z  ))\n",
    "b"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 81,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}\\frac{59 z}{34}\\\\\\frac{33 z}{34}\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "⎡59⋅z⎤\n",
       "⎢────⎥\n",
       "⎢ 34 ⎥\n",
       "⎢    ⎥\n",
       "⎢33⋅z⎥\n",
       "⎢────⎥\n",
       "⎣ 34 ⎦"
      ]
     },
     "execution_count": 81,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x = A.inv()*b\n",
    "x"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 82,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}5.726\\\\3.203\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "⎡5.726⎤\n",
       "⎢     ⎥\n",
       "⎣3.203⎦"
      ]
     },
     "execution_count": 82,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x.subs({z:3.3}).evalf(4)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 83,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "sympy.matrices.dense.MutableDenseMatrix"
      ]
     },
     "execution_count": 83,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "type(x)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "An alternative method of solving the same problem is to construct the system as a matrix in augmented form; that is the form obtained by appending the columns of (in our example) *A* and $\\vec{b}$ together. The augmented matrix is[1]:\n",
    "\n",
    "$$(A|\\vec{b})=\\left(\\begin{array}{cc|c}\n",
    "3 & 7 & 12z\\\\\n",
    "4 & -2 & 5z\\end{array} \\right),$$\n",
    "\n",
    " and as before we construct this as a SymPy `Matrix` object, but in this case we pass it to the `solve_linear_system()` function:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 84,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}3 & 7 & 12 z\\\\4 & -2 & 5 z\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "⎡3  7   12⋅z⎤\n",
       "⎢           ⎥\n",
       "⎣4  -2  5⋅z ⎦"
      ]
     },
     "execution_count": 84,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from sympy import Matrix, symbols, solve_linear_system\n",
    "x, y, z = symbols('x,y,z')\n",
    "system = Matrix(([3, 7, 12*z],[4, -2, 5*z]))\n",
    "system"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 85,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left\\{ x : \\frac{59 z}{34}, \\  y : \\frac{33 z}{34}\\right\\}$"
      ],
      "text/plain": [
       "⎧   59⋅z     33⋅z⎫\n",
       "⎨x: ────, y: ────⎬\n",
       "⎩    34       34 ⎭"
      ]
     },
     "execution_count": 85,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sol = solve_linear_system(system,x,y)\n",
    "sol"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 86,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "dict"
      ]
     },
     "execution_count": 86,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "type(sol)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 87,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "x = 5.726\n",
      "y = 3.203\n"
     ]
    }
   ],
   "source": [
    "# NBVAL_IGNORE_OUTPUT\n",
    "for k in sol.keys():\n",
    "    print(k,'=',sol[k].subs({z:3.3}).evalf(4))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A third option is the `solve()` method, whose arguments include the individual symbolic equations, rather than any matrices. Like `dsolve()` (see [ODEs](#Ordinary-differential-equations)), `solve()` expects either expressions which it will assume equal to zero, or `Equality` objects, which we can conveniently create with `Eq()`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 88,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left\\{ x : \\frac{59 z}{34}, \\  y : \\frac{33 z}{34}\\right\\}$"
      ],
      "text/plain": [
       "⎧   59⋅z     33⋅z⎫\n",
       "⎨x: ────, y: ────⎬\n",
       "⎩    34       34 ⎭"
      ]
     },
     "execution_count": 88,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from sympy import symbols,solve,Eq\n",
    "x, y, z = symbols('x,y,z')\n",
    "solve((Eq(3*x+7*y,12*z), Eq(4*x-2*y,5*z)), x, y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 89,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left\\{ x : \\frac{59 z}{34}, \\  y : \\frac{33 z}{34}\\right\\}$"
      ],
      "text/plain": [
       "⎧   59⋅z     33⋅z⎫\n",
       "⎨x: ────, y: ────⎬\n",
       "⎩    34       34 ⎭"
      ]
     },
     "execution_count": 89,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "solve((3*x+7*y-12*z, 4*x-2*y-5*z), x, y)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For more information, see `help(solve)` and `help(solve_linear_system)`.\n",
    "\n",
    "### Non linear equations\n",
    "\n",
    "Let’s solve a simple equation such as\n",
    "$x = x^2$. There are two obvious solutions: *x* = 0 and *x* = 1. How can we ask Sympy to compute these for us?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 90,
   "metadata": {},
   "outputs": [],
   "source": [
    "import sympy\n",
    "x, y, z = sympy.symbols('x, y, z')        # create some symbols\n",
    "eq = x - x ** 2                           # define the equation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 91,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[ 0, \\  1\\right]$"
      ],
      "text/plain": [
       "[0, 1]"
      ]
     },
     "execution_count": 91,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sympy.solve(eq, x)                        # solve eq = 0"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The `solve()` function expects an expression that as meant to be solve so that it evaluates to zero. For our example, we rewrite\n",
    "\n",
    "*x* = *x*<sup>2</sup>\n",
    " as\n",
    "*x* − *x*<sup>2</sup> = 0\n",
    " and then pass this to the solve function.\n",
    "\n",
    "Let’s repeat the same for the equation:\n",
    "*x* = *x*<sup>3</sup>\n",
    " and solve"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 92,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[ -1, \\  0, \\  1\\right]$"
      ],
      "text/plain": [
       "[-1, 0, 1]"
      ]
     },
     "execution_count": 92,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "eq = x - x ** 3                           # define the equation\n",
    "sympy.solve(eq, x)                        # solve eq = 0"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Output: LaTeX interface and pretty-printing\n",
    "\n",
    "As is the case with many computer algebra systems, SymPy has the ability to format its output as LaTeX code, for easy inclusion into documents.\n",
    "\n",
    "At the start of this chapter, we called:\n",
    "\n",
    "```python\n",
    "sympy.init_printing()\n",
    "```\n",
    "\n",
    "Sympy detected that it was in Jupyter, and enabled Latex output. The Jupyter Notebook supports (some) Latex, so this gives us the nicely formatted output above.\n",
    "\n",
    "We can also see the plain text output from Sympy, and the raw Latex code it creates:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 93,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "y**2/x**3 - y/x**2 + 1/x + O(y**3)\n"
     ]
    }
   ],
   "source": [
    "print(series(1/(x+y), y, 0, 3))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 94,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\\frac{y^{2}}{x^{3}} - \\frac{y}{x^{2}} + \\frac{1}{x} + O\\left(y^{3}\\right)\n"
     ]
    }
   ],
   "source": [
    "print(latex(series(1/(x+y), y, 0, 3)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 95,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "$\\frac{y^{2}}{x^{3}} - \\frac{y}{x^{2}} + 1 / x + O\\left(y^{3}\\right)$\n"
     ]
    }
   ],
   "source": [
    "print(latex(series(1/(x+y), y, 0, 3), mode='inline'))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Be aware that in its default mode, `latex()` outputs code that requires the `amsmath` package to be loaded via a `\\backslashusepackage{amsmath}` command in the document preamble.\n",
    "\n",
    "SymPy also supports a “pretty print” (`pprint()`) output routine, which produces better-formatted text output than the default printing routine, as illustrated below. Note features such as the subscripts for array elements whose names are of the form `T_n`, the italicised constant *e*, vertically-centred dots for multiplication, and the nicely-formed matrix borders and fractions.\n",
    "\n",
    "<img src=\"static/images/sympy.png\" alt=\"Nicely-formatted output from pprint().\" style=\"width:60.0%\" />\n",
    "\n",
    "Finally, SymPy offers `preview()`, which displays rendered output on screen (check `help(preview)` for details)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Automatic generation of C code\n",
    "\n",
    "A strong point of many symbolic libraries is that they can convert the symbolic expressions to C-code (or other code) that can subsequently be compiled for high execution speed. Here is an example that demonstrates this:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 96,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle x - \\frac{x^{3}}{6} + \\frac{x^{5}}{120} + O\\left(x^{6}\\right)$"
      ],
      "text/plain": [
       "     3     5        \n",
       "    x     x     ⎛ 6⎞\n",
       "x - ── + ─── + O⎝x ⎠\n",
       "    6    120        "
      ]
     },
     "execution_count": 96,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from sympy import *                                                                                    \n",
    "from sympy.utilities.codegen import codegen                                                                            \n",
    "x = Symbol('x')                                                                                                          \n",
    "sin(x).series(x, 0, 6)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 97,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "/******************************************************************************\n",
      " *                      Code generated with SymPy 1.11.1                      *\n",
      " *                                                                            *\n",
      " *              See http://www.sympy.org/ for more information.               *\n",
      " *                                                                            *\n",
      " *                       This file is part of 'project'                       *\n",
      " ******************************************************************************/\n",
      "#include \"taylor_sine.h\"\n",
      "#include <math.h>\n",
      "\n",
      "double taylor_sine(double x) {\n",
      "\n",
      "   double taylor_sine_result;\n",
      "   taylor_sine_result = x - 1.0/6.0*pow(x, 3) + (1.0/120.0)*pow(x, 5) + O(x**6);\n",
      "   return taylor_sine_result;\n",
      "\n",
      "}\n",
      "\n"
     ]
    }
   ],
   "source": [
    "# NBVAL_IGNORE_OUTPUT\n",
    "print(codegen((\"taylor_sine\",sin(x).series(x,0,6)), language='C')[0][1])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Related tools\n",
    "\n",
    "It is worth noting that the SAGE initiative <https://www.sagemath.org/> is trying to “create a viable free open source alternative to Magma, Maple, Mathematica and Matlab.” and includes the SymPy library among many others. Its symbolic capabilities are more powerful than SymPy’s, and SAGE, but the SymPy features will already cover many of the needs arising in science and engineering. SAGE includes the computer algebra system Maxima, which is also available standalone from <https://doc.sagemath.org/html/en/reference/interfaces/sage/interfaces/maxima_abstract.html>."
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "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.9.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}