Simpson's rule even number of points

Simpson_rule_even_points.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "f6cf89b9",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from scipy.integrate import simpson, quadrature"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "8fbe3a9f",
   "metadata": {},
   "outputs": [],
   "source": [
    "# equation from https://en.wikipedia.org/wiki/Simpson%27s_rule#Composite_Simpson's_rule_for_irregularly_spaced_data\n",
    "def remainder(y, x):\n",
    "    assert len(x) == 3\n",
    "    h = np.diff(x)\n",
    "    assert len(h) == 2\n",
    "    \n",
    "    alpha = 2*h[1]**2 + 3*h[0]*h[1]\n",
    "    alpha /= 6 * (h[1] + h[0])\n",
    "    \n",
    "    beta = h[1]**2 + 3 * h[0] * h[1]\n",
    "    beta /= 6 * h[0]\n",
    "    \n",
    "    nu = h[1]**3\n",
    "    nu /= 6 * h[0] * (h[0] + h[1])\n",
    "    \n",
    "    return alpha*y[-1] + beta*y[-2] - nu*y[-3]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "235418e9",
   "metadata": {},
   "source": [
    "## test out on quadratic, which should be handled exactly by Simpsons rule\n",
    "Definite integral should be 21"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "78844fc6",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "21.166666666666664"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "\n",
    "x = np.linspace(1, 4, 4)\n",
    "f = lambda x: x**2\n",
    "\n",
    "# this shows existing behaviour isn't exact.\n",
    "sim = simpson(f(x), x=x, even='avg')\n",
    "sim"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "03f7a480",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "21.0"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# this cell shows the integration is exact for an even number of points\n",
    "s = simpson(f(x[:-1]), x=x[:-1])\n",
    "r = remainder(f(x[-3:]), x=x[-3:])\n",
    "\n",
    "new = s+r\n",
    "new"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "e9a6d91b",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(21.0, 3.552713678800501e-15)"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ex = quadrature(f, 1, 4)\n",
    "ex"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "24d66b54",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(0.0, 0.7936507936507824)"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Percentage differences for new and existing behaviour\n",
    "100 * (new - ex[0])/ex[0], 100 * (sim - ex[0])/ex[0]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b9629ec4",
   "metadata": {},
   "source": [
    "## test out on np.cos"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "c9753738",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.3617683041940642"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# See what happens with a trig function.\n",
    "# first two terms of cos expansion is quadratic.\n",
    "x = np.linspace(0.5, 1, 4)\n",
    "f = lambda x: np.cos(x)\n",
    "\n",
    "sim = simpson(f(x), x=x, even='avg')\n",
    "sim"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "b564a2a3",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.3620708002437929"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# the improved sum\n",
    "s = simpson(f(x[:-1]), x=x[:-1])\n",
    "r = remainder(f(x[-3:]), x=x[-3:])\n",
    "\n",
    "new = s+r\n",
    "new"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "e1a9f644",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(0.36204544620289103, 2.8308679289601457e-09)"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ex = quadrature(f, 0.5, 1)\n",
    "ex"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "168d1e15",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(0.007002999531628406, -0.07654895586548137)"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Percentage differences for new and existing behaviour\n",
    "100 * (new - ex[0])/ex[0], 100 * (sim - ex[0])/ex[0]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "04d61dd0",
   "metadata": {},
   "source": [
    "## test out on np.sin"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "245d5c14",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.3370220715071605"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# See what happens with sin function.\n",
    "x = np.linspace(0.5, 1, 4)\n",
    "f = lambda x: np.sin(x)\n",
    "\n",
    "sim = simpson(f(x), x=x, even='avg')\n",
    "sim"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "id": "baed5ed9",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.33726008837298327"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# the improved sum\n",
    "s = simpson(f(x[:-1]), x=x[:-1])\n",
    "r = remainder(f(x[-3:]), x=x[-3:])\n",
    "\n",
    "new = s+r\n",
    "new"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "id": "0f4d2c25",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(0.3372802560214855, 2.637226437229856e-09)"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ex = quadrature(f, 0.5, 1)\n",
    "ex"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "id": "6d549bbd",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(-0.005979492763705395, -0.07654895586551003)"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Percentage differences for new and existing behaviour\n",
    "100 * (new - ex[0])/ex[0], 100 * (sim - ex[0])/ex[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "64f2af91",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.9.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}