{ "cells": [ { "cell_type": "markdown", "id": "6b37dc08", "metadata": { "id": "Z90IOhfZWL-Z" }, "source": [ " \n", "\n", "# **Numerical integration**\n", "\n", "An integral can be interpreted as the \"area under the curve\" of a function, or a generalization of area. Integrals have many applications in engineering and physics. However, it's sometimes very difficult to find the integral of a function in closed form (in terms of known functions). Therefore, in practice, we use various techniques for numerical integration to approximate this values." ] }, { "cell_type": "markdown", "id": "3eb71504", "metadata": { "id": "a0y5z-LCKGHc" }, "source": [ "## Integration formulas\n", "\n", "To find the numerical approximation of the integral of the function $f(x)$ in the interval [a,b], we should first discretize the interval into $n$ smaller subintervals of length $h=\\frac{b-a}{n}$. In each subinterval the integral is then computed using an estimate for the function value there. The main three approaches for finding the estimate of the function in each subinterval are (also shown in the figure below):\n", "\n", "\n", "1. **Midpoint/Rectangle rule**: the integral value in each interval is approximated by the area of the rectangle in each subinterval with width $h$ and a height equal to the value of the function in the midpoint of the subinterval: $\\int_{a}^{b}f(x)dx=\\sum_{i=0}^{n-1}hf(\\frac{x_i+x_{i+1}}{2})$\n", "\n", "2. **Trapezoidal rule**: the integral value in each interval is approximated by the area of the trapezoid in each subinterval with corners at the begining and the end of the interval: \n", "$\\int_{a}^{b}f(x)dx=\\sum_{i=0}^{n-1}h\\frac{f(x_i)+f(x_{i+1})}{2}$.\n", "\n", "3. **Simpson's rule**: Simpson’s rule approximates the area under two subintervals by fitting a second-order polynomial through the points and calculating its integral.\n", "\n", "![Integration formulas](../../images/integration-formulas.png)\n", "

" ] }, { "cell_type": "markdown", "id": "9e987e09", "metadata": { "id": "Xpb6CYWWUSWo" }, "source": [ "### Example 1: Using the Trapezoid rule\n", "The code below uses the Trapezoid rule to find $\\int_{0}^{2} \\cos (2\\pi x)dx$ with 21 evenly spaced grid points over the whole integration interval. How does it compare to the exact value?" ] }, { "cell_type": "code", "execution_count": 1, "id": "2cd8d20e", "metadata": { "colab": { "base_uri": "https://localhost:8080/" }, "id": "8rfay8mJUegW", "outputId": "65c67896-4132-4dcf-b2a8-edaf340165c3" }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "-6.66133814775094e-17\n" ] } ], "source": [ "import numpy as np\n", "\n", "a = 0\n", "b = 2\n", "n = 21 -1\n", "h = (b - a) / n\n", "x = np.linspace(a, b, n+1)\n", "f = np.cos(2*np.pi*x)\n", "\n", "I_trap = h * 0.5* np.sum ((f[0:-1]+f[1:]))\n", "print(I_trap)" ] }, { "cell_type": "markdown", "id": "4d7daeef", "metadata": { "id": "hOq3tVfFXhNC" }, "source": [ "## Numerical integration with SciPy\n", "\n", "You can also use the function [`scipy.integrate.quad(func, a, b)`](https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.quad.html) to find the numerical approximation of the integral of the function `func` in the interval [a,b].\n", "\n", "Below is an example to find $\\int_{0}^{2} \\cos (2\\pi x)dx$." ] }, { "cell_type": "code", "execution_count": 2, "id": "aee5ed0b", "metadata": { "colab": { "base_uri": "https://localhost:8080/" }, "id": "aEupVjb-U0dx", "outputId": "7fad99da-2add-4995-a052-b340cce63358" }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "The numerical result is 0.000000 (+-1.46797e-09)\n" ] } ], "source": [ "import numpy as np\n", "from scipy.integrate import quad\n", "\n", "# defining the function to integrate\n", "def f(x):\n", " return np.cos(2* np.pi * x )\n", "\n", "# finding integral using quad\n", "res, err = quad(f, 0, 2)\n", "\n", "print(\"The numerical result is {:f} (+-{:g})\"\n", " .format(res, err))" ] }, { "cell_type": "markdown", "id": "45fd4b55", "metadata": { "id": "rUIqmJfm5hp3" }, "source": [ "### Example 2: \n", "This assignment is dapted from Problem 6.41 in _Fundamentals of Electric Circuits_ (Alexander & Sadiku, 2017). \n", "\n", " The voltage across a $2H$ inductor is $20(1-\\exp^{-2t}) V$. If the initial current through the inductor is 0.3 A, find the current through the inductor at $t=1$s.\n", "Hint: $i(t)=\\frac{1}{L}\\int_{0}^{t}v(t) dt+i(0)$." ] }, { "cell_type": "code", "execution_count": 3, "id": "cdac0026", "metadata": { "colab": { "base_uri": "https://localhost:8080/" }, "id": "PAEDZqNCycHE", "outputId": "ec9b6791-4879-4f0c-8f95-11076c0ecac2" }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "5.9766764161830634\n" ] } ], "source": [ "import numpy as np\n", "from scipy.integrate import quad\n", "\n", "L= 2\n", "i0=0.3\n", "\n", "# defining the function to integrate\n", "def v(t):\n", " return 20*(1-np.exp(-2*t))\n", "\n", "# finding integral using quad\n", "res, err = quad(v, 0,1)\n", "i_1 = (1/L)*res+ i0\n", "\n", "print(i_1)" ] }, { "cell_type": "code", "execution_count": null, "id": "30feca00", "metadata": { "id": "TH-UpL27PvDZ" }, "outputs": [], "source": [] } ], "metadata": { "jupytext": { "text_representation": { "extension": ".md", "format_name": "myst", "format_version": 0.13, "jupytext_version": "1.16.4" } }, "kernelspec": { "display_name": "Python 3", "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.12.11" }, "source_map": [ 13, 21, 38, 43, 63, 71, 92, 100, 126 ] }, "nbformat": 4, "nbformat_minor": 5 }