{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "d7a8a988",
   "metadata": {},
   "source": [
    "# Sympy, Scipy, and integration"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "72789b15",
   "metadata": {},
   "source": [
    "The purpose of this notebook is to illustrate the application of sympy to make a numerical quadrature rule."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "239e6c81",
   "metadata": {},
   "source": [
    "Quadrature rules are applied to numerically approximate definite integrals."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8a39a619",
   "metadata": {},
   "source": [
    "The word *integration* applies to the mathematical problem, but in this course on software, mostly to the integration of several software packages."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "eca9fee3",
   "metadata": {},
   "source": [
    "### 1. Problem Statement"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8ce77baf",
   "metadata": {},
   "source": [
    "We show how to derivate an integration rule with sympy."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b0659e75",
   "metadata": {},
   "source": [
    "Our rule approximates an integral over [a,b]:\n",
    "\n",
    "$$\n",
    "    \\int_a^b f(x) dx = w_a f(a) + w_m f\\left( \\frac{a+b}{2} \\right) + w_b f(b)\n",
    "$$\n",
    "    \n",
    "evaluating a function $f$ at $a$, $m = (a+b)/2$, and $b$,\n",
    "so we look for 3 weights: $w_a$, $w_m$, and $w_b$."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "37a349f9",
   "metadata": {},
   "source": [
    "With 3 unknown weights we can require that all quadrics\n",
    "are integrated correctly because:\n",
    "    \n",
    "1. the integration operator is a linear operator and\n",
    "\n",
    "2. it suffices to require that the 3 basic functions 1, $x$, and $x^2$ are integrated correctly."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "27af1a49",
   "metadata": {},
   "source": [
    "### 2. A symbolic rule"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3fda4ab3",
   "metadata": {},
   "source": [
    "In Symbolic Computing we compute with symbols.  With ``sympy`` we declare the variables with ``var`` and define the quadrature rule for a symbolic function using ``Function``."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "19928220",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "making the rule WA*F(A) + WB*F(B) + WM*F(A/2 + B/2)\n"
     ]
    }
   ],
   "source": [
    "from sympy import var, Function\n",
    "A, B, WA, WM, WB = var('A, B, WA, WM, WB')\n",
    "RULE = lambda f: WA*f(A) + WM*f((A+B)/2) + WB*f(B)\n",
    "F = Function('F')\n",
    "print('making the rule', RULE(F))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2836c4ea",
   "metadata": {},
   "source": [
    "The requirement that the basic functions 1, $x$, and $x^2$ are integrated correctly leads to equations, which are the conditions on the three unknown weights $w_a$, $w_m$, $w_b$, respectively represented by ``WA``, ``WM``, and ``WB``."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "92f0e06f",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "solving 3 equations :\n",
      "A - B + WA + WB + WM == 0\n",
      "A**2/2 + A*WA - B**2/2 + B*WB + WM*(A/2 + B/2) == 0\n",
      "A**3/3 + A**2*WA - B**3/3 + B**2*WB + WM*(A/2 + B/2)**2 == 0\n"
     ]
    }
   ],
   "source": [
    "B0 = lambda x: 1\n",
    "B1 = lambda x: x\n",
    "B2 = lambda x: x**2\n",
    "# require that B0, B1, B2 are integrated exactly\n",
    "from sympy import integrate\n",
    "X = var('X')\n",
    "E0 = RULE(B0) - integrate(B0(X), (X, A, B))\n",
    "E1 = RULE(B1) - integrate(B1(X), (X, A, B))\n",
    "E2 = RULE(B2) - integrate(B2(X), (X, A, B))\n",
    "# we print the equations\n",
    "print('solving 3 equations :')\n",
    "print(E0, '== 0')\n",
    "print(E1, '== 0')\n",
    "print(E2, '== 0')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1bf7dfd6",
   "metadata": {},
   "source": [
    "In the cell above, observe that the basis functions are defined for the formal variable ``x``.  The functions are turned into expressions by evaluation at the letter ``X``.  The letter ``X`` is a global object to the entire session.  The formal variable ``x`` is local and makes sense only within the definition of the basis function.  Evaluaton of a function in a symbolic variable turns the function into an expression, suitable for input to ``integrate``."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5219450b",
   "metadata": {},
   "source": [
    "In solving the equations, we specify which are the variables of interest.  In this example, the left and right bounds of the integration interval are parameters.  We solve for the weights of the rule."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "d722b056",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "{WA: -A/6 + B/6, WM: -2*A/3 + 2*B/3, WB: -A/6 + B/6}\n"
     ]
    }
   ],
   "source": [
    "from sympy import solve, Subs, factor\n",
    "R = solve((E0, E1, E2), (WA, WM, WB))\n",
    "print(R)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "97edb06b",
   "metadata": {},
   "source": [
    "The ``solve`` returns a dictionary, assigned to ``R``, which has as its keys the variables given to ``solve`` and as values the computed solutions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "9ab8d78b",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Simpson formula : -(A - B)*(F(A) + F(B) + 4*F(A/2 + B/2))/6\n"
     ]
    }
   ],
   "source": [
    "V = RULE(F)\n",
    "S = Subs(V, (WA, WM, WB), (R[WA], R[WM], R[WB])).doit()\n",
    "FORMULA = factor(S)\n",
    "print('Simpson formula :', FORMULA)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2d214ca9",
   "metadata": {},
   "source": [
    "We computed the quadrature rule known as Simpson's formula:\n",
    "\n",
    "$$\n",
    "    \\int_a^b f(x) dx \\approx \\frac{b-a}{6} \\left( f(a) + 4 f\\left( \\frac{a+b}{2} \\right) + f(b) \\right).\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0837ce1a",
   "metadata": {},
   "source": [
    "### 3. Numerical Application"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "48f7d4d0",
   "metadata": {},
   "source": [
    "The rule, formally defined in ``FORMULA`` can be turned into an actual function, ready to approximate numerical integrals."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "496af76f",
   "metadata": {},
   "outputs": [],
   "source": [
    "from sympy import lambdify\n",
    "SIMPSON = lambdify((F, A, B), FORMULA)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "96ad7652",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "6.333333333333333"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "SIMPSON(lambda x: x**2, 2, 3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "d824336e",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\frac{19}{3}$"
      ],
      "text/plain": [
       "19/3"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "integrate(X**2, (X, 2, 3))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "bb47ca70",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "6.333333333333333"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "19/3"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c707e025",
   "metadata": {},
   "source": [
    "Let us verify if every quadric is integrated exactly."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "5b70cf44",
   "metadata": {},
   "outputs": [],
   "source": [
    "C0, C1, C2 = var('C0, C1, C2')\n",
    "QUADRIC = lambda x: C0 + C1*x + C2*x**2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "c164e8bc",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "approximation : -0.333333333333333*A**3*C2 - 0.5*A**2*C1 - 1.0*A*C0 + 0.333333333333333*B**3*C2 + 0.5*B**2*C1 + 1.0*B*C0\n",
      "  exact value : -A**3*C2/3 - A**2*C1/2 - A*C0 + B**3*C2/3 + B**2*C1/2 + B*C0\n",
      "    the error : 0\n"
     ]
    }
   ],
   "source": [
    "from sympy import simplify\n",
    "APPROX = simplify(SIMPSON(QUADRIC, A, B))\n",
    "print('approximation :', APPROX)\n",
    "EXACT = integrate(QUADRIC(X), (X, A, B))\n",
    "print('  exact value :', EXACT)\n",
    "print('    the error :', APPROX - EXACT)"
   ]
  }
 ],
 "metadata": {
  "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.10.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}