{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "cultural-moses",
   "metadata": {
    "jupyter": {
     "source_hidden": true
    }
   },
   "source": [
    "# Demo of Python in Jupyter Notebook"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "driving-table",
   "metadata": {},
   "outputs": [],
   "source": [
    "%%html\n",
    "<iframe width=\"560\" height=\"315\" src=\"https://www.youtube.com/embed/HW29067qVWk\" frameborder=\"0\" allow=\"autoplay; encrypted-media\" allowfullscreen></iframe>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "weird-fleece",
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "N = 50\n",
    "x = np.random.rand(N)\n",
    "y = np.random.rand(N)\n",
    "colors = np.random.rand(N)\n",
    "area = np.pi * (15 * np.random.rand(N))**2\n",
    "\n",
    "plt.scatter(x, y, s=area, c=colors, alpha=0.5)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "distinct-spell",
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import numpy as np\n",
    "\n",
    "df = pd.DataFrame(np.random.randn(10,5))\n",
    "df"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "answering-exhaust",
   "metadata": {},
   "source": [
    "## SymPy: Open Source Symbolic Mathematics\n",
    "\n",
    "This section uses the [SymPy](http://sympy.org) package to perform symbolic manipulations,\n",
    "and combined with numpy and matplotlib, also displays numerical visualizations of symbolically\n",
    "constructed expressions.\n",
    "\n",
    "We first load sympy printing extensions, as well as all of sympy:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "incorrect-royalty",
   "metadata": {},
   "outputs": [],
   "source": [
    "from __future__ import division\n",
    "from IPython.display import display\n",
    "\n",
    "from sympy.interactive import printing\n",
    "printing.init_printing(use_latex='mathjax')\n",
    "\n",
    "import sympy as sym\n",
    "from sympy import *\n",
    "x, y, z = symbols(\"x y z\")\n",
    "k, m, n = symbols(\"k m n\", integer=True)\n",
    "f, g, h = map(Function, 'fgh')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "consecutive-transsexual",
   "metadata": {},
   "source": [
    "### Elementary operations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "careful-minneapolis",
   "metadata": {},
   "outputs": [],
   "source": [
    "Rational(3,2)*pi + exp(I*x) / (x**2 + y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "intermediate-forest",
   "metadata": {},
   "outputs": [],
   "source": [
    "exp(I*x).subs(x,pi).evalf()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "continued-blackberry",
   "metadata": {},
   "outputs": [],
   "source": [
    "e = x + 2*y"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "advised-supplier",
   "metadata": {},
   "outputs": [],
   "source": [
    "srepr(e)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cordless-beginning",
   "metadata": {},
   "outputs": [],
   "source": [
    "exp(pi * sqrt(163)).evalf(50)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "conscious-neighborhood",
   "metadata": {},
   "source": [
    "### Algebra"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "impressive-white",
   "metadata": {},
   "outputs": [],
   "source": [
    "eq = ((x+y)**2 * (x+1))\n",
    "eq"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "electoral-vancouver",
   "metadata": {},
   "outputs": [],
   "source": [
    "expand(eq)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "literary-inflation",
   "metadata": {},
   "outputs": [],
   "source": [
    "a = 1/x + (x*sin(x) - 1)/x\n",
    "a"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "persistent-pledge",
   "metadata": {},
   "outputs": [],
   "source": [
    "simplify(a)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "configured-prerequisite",
   "metadata": {},
   "outputs": [],
   "source": [
    "eq = Eq(x**3 + 2*x**2 + 4*x + 8, 0)\n",
    "eq"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "rental-farming",
   "metadata": {},
   "outputs": [],
   "source": [
    "solve(eq, x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "minimal-archive",
   "metadata": {},
   "outputs": [],
   "source": [
    "a, b = symbols('a b')\n",
    "Sum(6*n**2 + 2**n, (n, a, b))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "exceptional-cathedral",
   "metadata": {},
   "source": [
    "### Calculus"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "pregnant-relaxation",
   "metadata": {},
   "outputs": [],
   "source": [
    "limit((sin(x)-x)/x**3, x, 0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "soviet-juice",
   "metadata": {},
   "outputs": [],
   "source": [
    "(1/cos(x)).series(x, 0, 6)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "breathing-intent",
   "metadata": {},
   "outputs": [],
   "source": [
    "diff(cos(x**2)**2 / (1+x), x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "noticed-drink",
   "metadata": {},
   "outputs": [],
   "source": [
    "integrate(x**2 * cos(x), (x, 0, pi/2))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "driven-corrections",
   "metadata": {},
   "outputs": [],
   "source": [
    "eqn = Eq(Derivative(f(x),x,x) + 9*f(x), 1)\n",
    "display(eqn)\n",
    "dsolve(eqn, f(x))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "recreational-manor",
   "metadata": {},
   "source": [
    "## Illustrating Taylor series\n",
    "\n",
    "We will define a function to compute the Taylor series expansions of a symbolically defined expression at\n",
    "various orders and visualize all the approximations together with the original function"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "visible-sociology",
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "interesting-equivalent",
   "metadata": {},
   "outputs": [],
   "source": [
    "# You can change the default figure size to be a bit larger if you want,\n",
    "# uncomment the next line for that:\n",
    "#plt.rc('figure', figsize=(10, 6))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "alert-treasury",
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_taylor_approximations(func, x0=None, orders=(2, 4), xrange=(0,1), yrange=None, npts=200):\n",
    "    \"\"\"Plot the Taylor series approximations to a function at various orders.\n",
    "\n",
    "    Parameters\n",
    "    ----------\n",
    "    func : a sympy function\n",
    "    x0 : float\n",
    "      Origin of the Taylor series expansion.  If not given, x0=xrange[0].\n",
    "    orders : list\n",
    "      List of integers with the orders of Taylor series to show.  Default is (2, 4).\n",
    "    xrange : 2-tuple or array.\n",
    "      Either an (xmin, xmax) tuple indicating the x range for the plot (default is (0, 1)),\n",
    "      or the actual array of values to use.\n",
    "    yrange : 2-tuple\n",
    "      (ymin, ymax) tuple indicating the y range for the plot.  If not given,\n",
    "      the full range of values will be automatically used. \n",
    "    npts : int\n",
    "      Number of points to sample the x range with.  Default is 200.\n",
    "    \"\"\"\n",
    "    if not callable(func):\n",
    "        raise ValueError('func must be callable')\n",
    "    if isinstance(xrange, (list, tuple)):\n",
    "        x = np.linspace(float(xrange[0]), float(xrange[1]), npts)\n",
    "    else:\n",
    "        x = xrange\n",
    "    if x0 is None: x0 = x[0]\n",
    "    xs = sym.Symbol('x')\n",
    "    # Make a numpy-callable form of the original function for plotting\n",
    "    fx = func(xs)\n",
    "    f = sym.lambdify(xs, fx, modules=['numpy'])\n",
    "    # We could use latex(fx) instead of str(), but matploblib gets confused\n",
    "    # with some of the (valid) latex constructs sympy emits.  So we play it safe.\n",
    "    plt.plot(x, f(x), label=str(fx), lw=2)\n",
    "    # Build the Taylor approximations, plotting as we go\n",
    "    apps = {}\n",
    "    for order in orders:\n",
    "        app = fx.series(xs, x0, n=order).removeO()\n",
    "        apps[order] = app\n",
    "        # Must be careful here: if the approximation is a constant, we can't\n",
    "        # blindly use lambdify as it won't do the right thing.  In that case, \n",
    "        # evaluate the number as a float and fill the y array with that value.\n",
    "        if isinstance(app, sym.core.numbers.Number):\n",
    "            y = np.zeros_like(x)\n",
    "            y.fill(app.evalf())\n",
    "        else:\n",
    "            fa = sym.lambdify(xs, app, modules=['numpy'])\n",
    "            y = fa(x)\n",
    "        tex = sym.latex(app).replace('$', '')\n",
    "        plt.plot(x, y, label=r'$n=%s:\\, %s$' % (order, tex) )\n",
    "        \n",
    "    # Plot refinements\n",
    "    if yrange is not None:\n",
    "        plt.ylim(*yrange)\n",
    "    plt.grid()\n",
    "    plt.legend(loc='best').get_frame().set_alpha(0.8)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "surprising-start",
   "metadata": {},
   "source": [
    "With this function defined, we can now use it for any sympy function or expression"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "spoken-exhibition",
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_taylor_approximations(sin, 0, [2, 4, 6], (0, 2*pi), (-2,2))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "behind-scout",
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_taylor_approximations(cos, 0, [2, 4, 6], (0, 2*pi), (-2,2))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "floral-paint",
   "metadata": {},
   "source": [
    "This shows easily how a Taylor series is useless beyond its convergence radius, illustrated by a simple function that has singularities on the real axis:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "aerial-tennis",
   "metadata": {},
   "outputs": [],
   "source": [
    "# For an expression made from elementary functions, we must first make it into\n",
    "# a callable function, the simplest way is to use the Python lambda construct.\n",
    "plot_taylor_approximations(lambda x: 1/cos(x), 0, [2,4,6], (0, 2*pi), (-5,5))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "raising-raleigh",
   "metadata": {},
   "source": [
    "## Basic Numerical Integration: the Trapezoid Rule"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "wound-perspective",
   "metadata": {},
   "source": [
    "A simple illustration of the trapezoid rule for definite integration:\n",
    "\n",
    "$$\n",
    "\\int_{a}^{b} f(x)\\, dx \\approx \\frac{1}{2} \\sum_{k=1}^{N} \\left( x_{k} - x_{k-1} \\right) \\left( f(x_{k}) + f(x_{k-1}) \\right).\n",
    "$$\n",
    "<br>\n",
    "First, we define a simple function and sample it between 0 and 10 at 200 points"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "brown-denmark",
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "requested-convert",
   "metadata": {},
   "outputs": [],
   "source": [
    "def f(x):\n",
    "    return (x-3)*(x-5)*(x-7)+85\n",
    "\n",
    "x = np.linspace(0, 10, 200)\n",
    "y = f(x)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "incorrect-passenger",
   "metadata": {},
   "source": [
    "Choose a region to integrate over and take only a few points in that region"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "iraqi-navigation",
   "metadata": {},
   "outputs": [],
   "source": [
    "a, b = 1, 8 # the left and right boundaries\n",
    "N = 5 # the number of points\n",
    "xint = np.linspace(a, b, N)\n",
    "yint = f(xint)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "considered-brick",
   "metadata": {},
   "source": [
    "Plot both the function and the area below it in the trapezoid approximation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "assumed-practitioner",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(x, y, lw=2)\n",
    "plt.axis([0, 9, 0, 140])\n",
    "plt.fill_between(xint, 0, yint, facecolor='gray', alpha=0.4)\n",
    "plt.text(0.5 * (a + b), 30,r\"$\\int_a^b f(x)dx$\", horizontalalignment='center', fontsize=20);"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "antique-thomson",
   "metadata": {},
   "source": [
    "Compute the integral both at high accuracy and with the trapezoid approximation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "patent-psychology",
   "metadata": {},
   "outputs": [],
   "source": [
    "from __future__ import print_function\n",
    "from scipy.integrate import quad\n",
    "integral, error = quad(f, a, b)\n",
    "integral_trapezoid = sum( (xint[1:] - xint[:-1]) * (yint[1:] + yint[:-1]) ) / 2\n",
    "print(\"The integral is:\", integral, \"+/-\", error)\n",
    "print(\"The trapezoid approximation with\", len(xint), \"points is:\", integral_trapezoid)"
   ]
  }
 ],
 "metadata": {
  "@webio": {
   "lastCommId": null,
   "lastKernelId": null
  },
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.13.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
