{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# SymPy -- symbolic computation\n",
    "\n",
    "Symbolic computation is the computation of mathematical objects in an analytical way. \n",
    "The mathematical objects are represented exactly and not only approximately. For example $\\sqrt{8}$ will not approximated as $\\sqrt{8} = 2.82842712475$, it will be done as:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle 2 \\sqrt{2}$"
      ],
      "text/plain": [
       "2*sqrt(2)"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import sympy as sp\n",
    "sp.sqrt(8)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Furthermore one can define expressions, equations and more "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "x + 3*y\n"
     ]
    },
    {
     "data": {
      "text/latex": [
       "$\\displaystyle x + 3 y + 1$"
      ],
      "text/plain": [
       "x + 3*y + 1"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x, y = sp.symbols('x y')\n",
    "expr=x+3*y\n",
    "print(expr)\n",
    "expr+1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Differentiation, Integration, computation of limits is possible:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle 6 x + \\cos{\\left(x \\right)} + 4$"
      ],
      "text/plain": [
       "6*x + cos(x) + 4"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sp.diff(sp.sin(x)+3*x**2+4*x,x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left(x^{2} - 2 x + 2\\right) e^{x}$"
      ],
      "text/plain": [
       "(x**2 - 2*x + 2)*exp(x)"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sp.integrate(x**2*sp.exp(x),x)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle -2 + 82 e^{10}$"
      ],
      "text/plain": [
       "-2 + 82*exp(10)"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sp.integrate(x**2*sp.exp(x),(x,0,10))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle 1806168.19517415$"
      ],
      "text/plain": [
       "1806168.19517415"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sp.integrate(x**2*sp.exp(x),(x,0,10)).evalf()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle 1$"
      ],
      "text/plain": [
       "1"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sp.limit(sp.sin(x)/x,x,0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Task\n",
    "\n",
    "Differentiate, integrate some expressions of your choice and compute the limits of some series"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "### Floating point representation\n",
    "\n",
    "Sometimes it is useful to let SymPy evaluate a floating point approximation (up to a user specified number of digits).\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle 1.414213562$"
      ],
      "text/plain": [
       "1.414213562"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sp.sqrt(2).evalf(10)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Task\n",
    "\n",
    "Evaluate $\\pi$ from sympy with hundred digits. \n",
    "\n",
    "Then write a functions which uses Leibniz Formula \n",
    "$$ 1 - \\frac{1}{3} + \\frac{1}{5} -  \\frac{1}{7}  + \\frac{1}{9}  + - \\cdots = \\frac{\\pi}{4} $$\n",
    "to approximate $\\pi$ with $n$ summands from the formula above. Make a convergence test - how many summands you need to get $m=5$ digits of $\\pi$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "# get 100 digits of pi "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0\n"
     ]
    }
   ],
   "source": [
    "def piapprox(n):\n",
    "    p=0\n",
    "    # TODO \n",
    "    return p\n",
    "\n",
    "print(4*piapprox(100000000))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Solving equations\n",
    "\n",
    "Sympy can solve equations, for example roots of polynomials"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[-1, 4]"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sp.solve(-x**2+3*x+4,x)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Not all is possible in exact computation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "ename": "NotImplementedError",
     "evalue": "multiple generators [exp(x), sin(x)]\nNo algorithms are implemented to solve equation exp(x) - sin(x)",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mNotImplementedError\u001b[0m                       Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-17-94031e630926>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[0;32m----> 1\u001b[0;31m \u001b[0msp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0msolve\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0msp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mexp\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mx\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m-\u001b[0m\u001b[0msp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0msin\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mx\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mx\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m",
      "\u001b[0;32m~/.local/lib/python3.6/site-packages/sympy/solvers/solvers.py\u001b[0m in \u001b[0;36msolve\u001b[0;34m(f, *symbols, **flags)\u001b[0m\n\u001b[1;32m   1172\u001b[0m     \u001b[0;31m###########################################################################\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   1173\u001b[0m     \u001b[0;32mif\u001b[0m \u001b[0mbare_f\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m-> 1174\u001b[0;31m         \u001b[0msolution\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0m_solve\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mf\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m*\u001b[0m\u001b[0msymbols\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m**\u001b[0m\u001b[0mflags\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m   1175\u001b[0m     \u001b[0;32melse\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   1176\u001b[0m         \u001b[0msolution\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0m_solve_system\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mf\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0msymbols\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m**\u001b[0m\u001b[0mflags\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/.local/lib/python3.6/site-packages/sympy/solvers/solvers.py\u001b[0m in \u001b[0;36m_solve\u001b[0;34m(f, *symbols, **flags)\u001b[0m\n\u001b[1;32m   1746\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   1747\u001b[0m     \u001b[0;32mif\u001b[0m \u001b[0mresult\u001b[0m \u001b[0;32mis\u001b[0m \u001b[0;32mFalse\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m-> 1748\u001b[0;31m         \u001b[0;32mraise\u001b[0m \u001b[0mNotImplementedError\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m'\\n'\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mjoin\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mmsg\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mnot_impl_msg\u001b[0m \u001b[0;34m%\u001b[0m \u001b[0mf\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m   1749\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   1750\u001b[0m     \u001b[0;32mif\u001b[0m \u001b[0mflags\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mget\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m'simplify'\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;32mTrue\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mNotImplementedError\u001b[0m: multiple generators [exp(x), sin(x)]\nNo algorithms are implemented to solve equation exp(x) - sin(x)"
     ]
    }
   ],
   "source": [
    "sp.solve(sp.exp(x)-sp.sin(x),x)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "but it can be approximated in a numerical way. Play around with different starting values."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle -15.7079634186507$"
      ],
      "text/plain": [
       "-15.7079634186507"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sp.nsolve(sp.exp(x)-sp.sin(x),0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It is possible to solve systems, for example the intersection of a circle and a line:\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[(1/5 - sqrt(179)/10, 1/10 + sqrt(179)/5),\n",
       " (1/5 + sqrt(179)/10, 1/10 - sqrt(179)/5)]"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "eq1 = x**2 + y**2 - 9  # circle of radius 3\n",
    "eq2 = 2*x + y - sp.Rational(1,2)   # line y(x) = -2*x + 1/2\n",
    "#eq2 = 2*x + y - 1/2  # try out what happens here\n",
    "\n",
    "sp.solve([eq1, eq2], [x, y])\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Task\n",
    "\n",
    "Make a geometrical interpretation to the following system of equations\n",
    "$$ x^2+y^2+z^2=1 $$\n",
    "$$ z=x^2+y^2 $$\n",
    "$$ y=0 $$\n",
    "and then solve this system, dump the solutions as expressions and as floating point numbers.\n",
    "Then interprete the solutions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [],
   "source": [
    "z=sp.symbols('z') # x and y are already symbols \n",
    "# TODO"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Task\n",
    "\n",
    "Now define your own geometrical problem in the space $R^3$. It should have $4$ points in $R^3$ as solution."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# TODO"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Task\n",
    "\n",
    "Now a task for integration and solving of equations. We consider the affine linear function\n",
    "$$y=m x + 1 \\quad m \\in R$$ \n",
    "and the parabola\n",
    "$$y=(x-1)^2$$\n",
    "\n",
    "Make a sketch of the functions (on a sheet of paper, not here) \n",
    "and then compute $m$ such that the area between the graphs of this functions is $1$.\n",
    "  * Hint 1 : you need the intersection points of the functions, depending on m\n",
    "  * Hint 2 : you get the area between the functions with integrals between the intersection points \n",
    "  * Hint 3 : then you have the area A(m) and can solve A(m)=1\n",
    "  * Hint 4 : check your solution at the end of all with a probe"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# TODO "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Differential equations \n",
    "\n",
    "Sympy can solve differential equations like\n",
    "\n",
    "$$u'(t)=a u(t) \\qquad a \\in R$$ \n",
    "\n",
    "which is an example for exponential growth (a>0) or exponential decay (a<0). Such equations can model \n",
    "the growth of microbes or radioactive disintegration."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle u{\\left(t \\right)} = C_{1} e^{a t}$"
      ],
      "text/plain": [
       "Eq(u(t), C1*exp(a*t))"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "u=sp.Function('u')\n",
    "t=sp.symbols('t')\n",
    "a=sp.symbols('a')\n",
    "sp.dsolve(u(t).diff(t, 1) - a*u(t))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Task \n",
    "\n",
    "Solve the problem for the logistic equation\n",
    "\n",
    "$$ u'(t) = au(t) - bu^2(t) $$\n",
    "\n",
    "and interpret the solution for $a=b=1$. Find out if there exist equilibrium states where $u(t)=const$ which is equivalent to $u'(t)=0$\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.6.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}