{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Compare numerial integral calculation\n",
    "==="
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "1. quadrature\n",
    "\n",
    "2. trapezoidal fomula\n",
    "\n",
    "3. simpson's rule"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def quad(func, lim=(0, 1), num=100, cof=1.0):\n",
    "    \"\"\"quadrature integral approximation\"\"\"\n",
    "    x = np.linspace(*lim, num + 1) * cof\n",
    "    k = np.ones(x.size) * (x[1] - x[0])\n",
    "    return k.dot(func(x))\n",
    "\n",
    "\n",
    "def trape(func, lim=(0, 1), num=100, cof=1.0):\n",
    "    \"\"\"trapezoidal integral approximation\"\"\"\n",
    "    x = np.linspace(*lim, num + 1) * cof\n",
    "    k = np.ones(x.size) * (x[1] - x[0])\n",
    "    k[0] *= 0.5\n",
    "    k[-1] *= 0.5\n",
    "    return k.dot(func(x))\n",
    "\n",
    "\n",
    "def simp(func, lim=(0, 1), num=100, cof=1.0):\n",
    "    \"\"\"simpson integral approximation\"\"\"\n",
    "    x = np.linspace(*lim, num + 1) * cof\n",
    "    k = np.ones(x.size) * (x[1] - x[0]) / 3.0\n",
    "    k[1:-1:2] *= 4.0\n",
    "    k[2:-2:2] *= 2.0\n",
    "    return k.dot(func(x))\n",
    "\n",
    "\n",
    "def sin_integ(lim=(0, 1), cof=1.0):\n",
    "    \"\"\"sin integral\"\"\"\n",
    "    return np.cos(cof * lim[0]) - np.cos(cof * lim[1])\n",
    "\n",
    "\n",
    "def error(ref, true):\n",
    "    \"\"\"relative error\"\"\"\n",
    "    return 100 * abs(ref - true) / true"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Numerical Integration of $f(x) = x^2$\n",
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The number of descritazation is 100 points along to x axis (0, 1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "true = 1.0 / 3.0  # analytical solution\n",
    "def func(x):\n",
    "    return x ** 2\n",
    "qua = quad(func)\n",
    "trap = trape(func)\n",
    "simpson = simp(func)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"quadrature relative error  : {:.2e} %\".format(error(qua, true)))\n",
    "print(\"trapezoidal relative error : {:.2e} %\".format(error(trap, true)))\n",
    "print(\"simpsons relative error    : {:.2e} %\".format(error(simpson, true)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Numerical Integration of $f(x) = \\sin x$\n",
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The number of descritazation is 100 points along to x, the range of which is (0, $\\pi$)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "true = 2\n",
    "def func(x):\n",
    "    return np.sin(x)\n",
    "qua = quad(func, lim=(0, np.pi))\n",
    "trap = trape(func, lim=(0, np.pi))\n",
    "simpson = simp(func, lim=(0, np.pi))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"quadrature relative error  : {:.2e} %\".format(error(qua, true)))\n",
    "print(\"trapezoidal relative error : {:.2e} %\".format(error(trap, true)))\n",
    "print(\"simpsons relative error    : {:.2e} %\".format(error(simpson, true)))"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "cherab-phix-dev",
   "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.15 | packaged by conda-forge | (main, Nov 22 2022, 15:55:03) \n[GCC 10.4.0]"
  },
  "vscode": {
   "interpreter": {
    "hash": "2725905a4c02db19e04df9b8fdbbe5ec65a73ea52bebaf9474aa1cc98819834c"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}