{ "cells": [ { "cell_type": "markdown", "id": "c78c3ac7-3476-4970-b33e-e1d3e6063257", "metadata": {}, "source": [ "# Gaussian quadrature" ] }, { "cell_type": "markdown", "id": "5916ba10-f7a4-477a-a2a1-32d26c4ac906", "metadata": {}, "source": [ "A quadrature rule is the approximation of a definite integral by a\n", "finite sum of the form\n", "\n", "$$\\int\\limits_{a}^{b} \\omega(x) \\, f(x) \\, \\mathrm{d}x \\approx \\sum_{i=1}^{n} w_i f(x_i) .$$\n", "\n", "Here $\\omega(x)$ is a chosen *weight function*, $x_1, \\ldots, x_n$ and $w_1, \\ldots, w_n$ are referred to as the\n", "quadrature nodes and weights, respectively. In 1814 Gauss described a\n", "choice for the nodes and weights that is optimal for $a = -1$, $b = 1$, $\\omega(x) = 1$ in the sense that for\n", "each $n$ it exactly integrates polynomials $f(x)$ up to degree $2n - 1$.\n", "\n", "\n", "We illustrate the idea of Gaussian quadrature using a simple examples." ] }, { "cell_type": "markdown", "id": "d635e193-b9ca-40a3-a447-c8c7c7853b02", "metadata": {}, "source": [ "## Gauss-Hermite quadrature" ] }, { "cell_type": "markdown", "id": "59924041-b326-47e3-9925-26687c4af855", "metadata": {}, "source": [ "Gauss-Hermite quadrature is a form of Gaussian quadrature for\n", "approximating the value of integrals of the following form:\n", "\n", "$$ \\int\\limits_{-\\infty}^{+\\infty} e^{-x^2} f(x) \\, \\mathrm{d}x .$$\n", "\n", "In this case \n", "\n", "$$\\int\\limits_{-\\infty}^{+\\infty} e^{-x^2} f(x) \\, \\mathrm{d}x \\approx \\sum _{i=1}^{n} w_i f(x_i) .$$\n", "\n", "where n is the number of nodes." ] }, { "cell_type": "markdown", "id": "902e4a0d-dc41-49b6-9d6a-25d765ffb523", "metadata": {}, "source": [ "Consider the two-point Gauss-Hermite quadrature\n", "\n", "$$\\int\\limits_{-\\infty}^{+\\infty} e^{-x^2} f(x) \\, \\mathrm{d}x \\approx w_1 f(x_1) + w_2 f(x_2) .$$ \n", "\n", "the integration range in Gauss-Hermite quadrature is symmetric with respect to the origin, so must be\n", "the coordinates of the nodes and the weights:\n", "\n", "$$x_1 = - x_2, \\qquad w_1 = w_2 .$$\n", "\n", "Thus two-point rule contains just two parameters:\n", "\n", "$$\\int\\limits_{-\\infty}^{+\\infty} e^{-x^2} f(x) \\, \\mathrm{d}x \\approx w f(x) + w f(-x) .$$\n", "\n", "Determine the unknown parameters from the requirements that\n", "the quadrature is exact for $f(x) = 1$ and $f(x) = x^2$:\n", "\n", "$$\\int\\limits_{-\\infty}^{+\\infty} e^{-x^2} \\, \\mathrm{d}x = 2w, \\qquad\n", " \\int\\limits_{-\\infty}^{+\\infty} e^{-x^2} x^2 \\, \\mathrm{d}x = 2 w x^2.$$\n", "\n", "Recall that \n", "$$\\int\\limits_{-\\infty}^{+\\infty} e^{-x^2} \\, \\mathrm{d}x = \\sqrt{\\pi}, \\qquad\n", "\\int\\limits_{-\\infty}^{+\\infty} e^{-x^2} x^2 = \\frac{\\sqrt{\\pi}}{2}.$$\n", "\n", "We obtain the following equations for $x$ and $w$:\n", "\n", "$$2w = \\sqrt{\\pi}, \\qquad 2 w x^2 = \\frac{\\sqrt{\\pi}}{2}$$\n", "\n", "Thus,\n", "\n", "$$w = \\frac{\\sqrt{\\pi}}{2}, \\qquad x = \\frac{1}{\\sqrt{2}},$$\n", "\n", "$$\\int\\limits_{-\\infty}^{+\\infty} e^{-x^2} f(x) \\, \\mathrm{d}x \\approx \\frac{\\sqrt{\\pi}}{2} \\left[\n", " f\\left(-\\frac{1}{\\sqrt{2}}\\right) + f\\left(\\frac{1}{\\sqrt{2}}\\right)\n", " \\right] .$$" ] }, { "cell_type": "markdown", "id": "dd965b6e-3b5c-45e8-a818-f12d174f5a6a", "metadata": {}, "source": [ "## Numerical examples" ] }, { "cell_type": "code", "execution_count": null, "id": "ec18db06-0703-4730-a2a3-c70b5c93ae64", "metadata": {}, "outputs": [], "source": [ "\n", "using QuadGK\n", "using PyPlot" ] }, { "cell_type": "markdown", "id": "755897cd-2d3b-4b9d-9118-490f404aa841", "metadata": {}, "source": [ "### Example 1" ] }, { "cell_type": "code", "execution_count": null, "id": "c9e68861-17c6-43c3-8d0e-29f15401af95", "metadata": {}, "outputs": [], "source": [ "\n", "f1(x) = cos(x/pi)" ] }, { "cell_type": "code", "execution_count": null, "id": "6d4bd98b-6c90-45bc-8204-e7b05e2c3d7f", "metadata": {}, "outputs": [], "source": [ "\n", "xx = range(-3.5, 3.5, 100)\n", "xs = range(-2.0, 2.0, 100);" ] }, { "cell_type": "code", "execution_count": null, "id": "4bed66a3-67aa-4db2-adb9-235b9901de1f", "metadata": {}, "outputs": [], "source": [ "\n", "plot(xs, f1.(xs), label=\"f1(x)\", color=\"black\")\n", "plot(xx, f1.(xx), linestyle=\"dotted\", color=\"black\")\n", "plot(xx, exp.(-xx .^ 2), label=L\"e^{-x^2}\", linestyle=\"dashed\")\n", "grid(true)\n", "legend();" ] }, { "cell_type": "markdown", "id": "a4d8c624-e635-4a78-8b14-cd46a5644e37", "metadata": {}, "source": [ "Two-point Gauss-Hermite quadrature:" ] }, { "cell_type": "code", "execution_count": null, "id": "46e4ec2c-65bd-4f07-954e-1b435503e8ed", "metadata": {}, "outputs": [], "source": [ "\n", "res1 = sqrt(pi) / 2 * (f1(1/sqrt(2)) + f1(-1/sqrt(2)))" ] }, { "cell_type": "markdown", "id": "7d49773f-17cb-478f-904f-468168491973", "metadata": {}, "source": [ "Adaptive integrator:" ] }, { "cell_type": "code", "execution_count": null, "id": "3e8b4969-58a1-49ce-9c50-816ae2d3b6a9", "metadata": {}, "outputs": [], "source": [ "\n", "g1(x) = f1(x) * exp(-x^2)\n", "chk1, error = quadgk(g1, -Inf, Inf)" ] }, { "cell_type": "markdown", "id": "7c5affa7-c7d7-4817-9947-c334c36dfa9f", "metadata": {}, "source": [ "Relative error:" ] }, { "cell_type": "code", "execution_count": null, "id": "c9c521e5-5df1-4aeb-8630-20deab74c762", "metadata": {}, "outputs": [], "source": [ "\n", "abs(res1 - chk1)/chk1" ] }, { "cell_type": "markdown", "id": "94024cc4-c507-48a3-845c-402dfca1a8aa", "metadata": {}, "source": [ "### Example 2" ] }, { "cell_type": "code", "execution_count": null, "id": "a10762da-9760-4c31-ba8b-2289e0d62f38", "metadata": {}, "outputs": [], "source": [ "\n", "f2(x) = 1/(1 + exp(x))" ] }, { "cell_type": "code", "execution_count": null, "id": "9c932e6d-5163-4c8d-ab6c-b43b623481de", "metadata": {}, "outputs": [], "source": [ "\n", "plot(xs, f2.(xs), label=\"f2(x)\", color=\"black\")\n", "plot(xx, f2.(xx), linestyle=\"dotted\", color=\"black\")\n", "plot(xx, exp.(-xx .^ 2), label=L\"e^{-x^2}\", linestyle=\"dashed\")\n", "grid(true)\n", "legend();" ] }, { "cell_type": "markdown", "id": "153fd6b2-72df-4d33-8783-8d86766ae05b", "metadata": {}, "source": [ "Two-point Gauss-Hermite quadrature:" ] }, { "cell_type": "code", "execution_count": null, "id": "4fdd5562-1bbf-4cdc-8533-343554fcf908", "metadata": {}, "outputs": [], "source": [ "\n", "res2 = sqrt(pi) / 2 * (f2(1/sqrt(2)) + f2(-1/sqrt(2)))" ] }, { "cell_type": "markdown", "id": "993bda86-3876-4b96-8641-e4164a0f6511", "metadata": {}, "source": [ "Adaptive integrator:" ] }, { "cell_type": "code", "execution_count": null, "id": "e3b7d6a6-3e26-4c46-9cde-6162b4ea0cab", "metadata": {}, "outputs": [], "source": [ "\n", "g2(x) = f2(x) * exp(-x^2)\n", "chk2, error = quadgk(g2, -Inf, Inf)" ] }, { "cell_type": "markdown", "id": "0bcfc8f8-d36f-4d8e-805d-56dae7229c3d", "metadata": {}, "source": [ "Relative error:" ] }, { "cell_type": "code", "execution_count": null, "id": "6ff4f28d-e488-44d0-a4e8-854e83f8d7ba", "metadata": {}, "outputs": [], "source": [ "\n", "abs(res2 - chk2)/chk2" ] }, { "cell_type": "code", "execution_count": null, "id": "86f9e27c-5952-4c57-8e66-5e38270953d2", "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "Julia 1.10.5", "language": "julia", "name": "julia-1.10" }, "language_info": { "file_extension": ".jl", "mimetype": "application/julia", "name": "julia", "version": "1.10.5" } }, "nbformat": 4, "nbformat_minor": 5 }