{ "cells": [ { "cell_type": "markdown", "id": "628f9256-8e52-4780-9d3e-b896e344b00c", "metadata": {}, "source": [ "# Trapezoidal integration formula and Richardson extrapolation" ] }, { "cell_type": "markdown", "id": "0c54ed57-45e1-4408-98d5-e0f1e7654873", "metadata": {}, "source": [ "
\n", "Trapezoidal rule with the integration step $h$:\n", "$$\n", "\\int_a^b f(x) \\, \\mathrm{d}x =\n", "\\frac{h}{2} \\left[\n", "f(x_1) + 2f(x_2) + 2f(x_3) + 2f(x_4)+ \\ldots + 2f(x_{n-1}) + f(x_n)\n", "\\right] + \\alpha h^2 + O(h^p) = T(h, n) + \\alpha h^2 + O(h^p) , \\tag{1}\n", "$$\n", "where for later use we defined\n", "$$T(h, n) \\equiv \\frac{h}{2} \\left[\n", "f(x_1) + 2f(x_2) + 2f(x_3) + 2f(x_4)+ \\ldots + 2f(x_{n-1}) + f(x_n)\n", "\\right] , \\tag{2}$$\n", "and\n", "$$\n", "h = \\frac{b-a}{n - 1},\n", "$$\n", "$$\n", "x_i = a + (i - 1) \\, h, \\quad i = 1, \\ldots, n,\n", "$$\n", "$\\alpha h^2$ is the leading discretization error term of the\n", "algorithm, the coefficient $\\alpha$ is unknown. The term $O(h^p)$\n", "indicates the higher order error terms, $p > 2$. \n", "
" ] }, { "cell_type": "code", "execution_count": null, "id": "dd5007fa-a3f9-45a0-b8e4-ea380c4ecb98", "metadata": {}, "outputs": [], "source": [ "using PyPlot" ] }, { "cell_type": "code", "execution_count": null, "id": "117777a7-8308-4286-bcef-85b40c127c50", "metadata": {}, "outputs": [], "source": [ "\"\"\"\n", " ans = mytrapezoids(fun, a, b, n)\n", "\n", "Numerically evaluate the integral int_a^b fun(x) dx using the\n", "trapezoidal rule: I = h/2*(f_1 + 2f_2 + ... + 2f_{n-1} + f_n),\n", "where h = (b - a)/(n - 1), x_i = a + (i - 1)*h, f_i = fun(x_i).\n", "\"\"\"\n", "function mytrapezoids(fun, a, b, n)\n", " h = (b - a)/(n - 1)\n", " s1 = fun(a) + fun(b)\n", " s2 = 0.0\n", " for i = 2:(n-1)\n", " s2 += fun(a + (i-1)*h)\n", " end\n", " return (s1 + 2*s2)*h/2\n", "end" ] }, { "cell_type": "markdown", "id": "d46c6b5e-38e8-43df-85a6-70a0ccf76b11", "metadata": {}, "source": [ "### Testing" ] }, { "cell_type": "code", "execution_count": null, "id": "23e405a1-484f-40f2-a5db-c689e552aeb4", "metadata": {}, "outputs": [], "source": [ "# Simple test 1\n", "fun1(x) = x\n", "a = -1.0\n", "b = 1.0\n", "np = 20\n", "exact = 0.0\n", "abs(mytrapezoids(fun1, a, b, np) - exact) " ] }, { "cell_type": "code", "execution_count": null, "id": "738e7463-2222-4a29-b60f-73dd15c4f1c2", "metadata": {}, "outputs": [], "source": [ "# Simple test 2\n", "fun2(x) = exp(x)\n", "a = 0.0\n", "b = 1.0\n", "np = 200\n", "exact = exp(1) - 1.0\n", "abs(mytrapezoids(fun2, a, b, np) - exact)" ] }, { "cell_type": "markdown", "id": "9b3f98df-c9ed-47fe-b9a9-a218da4d081f", "metadata": {}, "source": [ "### Error of trapezoidal formula" ] }, { "cell_type": "code", "execution_count": null, "id": "36816c8a-9754-4431-9dc6-88b806350ef3", "metadata": {}, "outputs": [], "source": [ "# An analysis of the error of trapezoidal rule vs integration step h\n", "fun0(x) = sin(x)\n", "a = 0.0\n", "b = pi\n", "exact = 2.0\n", "ndp = 10\n", "errt = zeros(ndp)\n", "hh = zeros(ndp);" ] }, { "cell_type": "code", "execution_count": null, "id": "2bb63042-e732-4adb-b0c8-079a3eddc792", "metadata": {}, "outputs": [], "source": [ "for k = 1:ndp \n", " np = 2^k + 1\n", " hh[k] = (b - a)/(np - 1)\n", " errt[k] = abs(mytrapezoids(fun0, a, b, np) - exact)\n", "end" ] }, { "cell_type": "code", "execution_count": null, "id": "737191c8-704f-4f31-8f35-bf0c8b306a59", "metadata": {}, "outputs": [], "source": [ "loglog(hh, errt, label=\"numerical experiment\", marker=\".\")\n", "loglog(hh, hh .^ 2, linestyle=\"dashed\", label=L\"hh^2\")\n", "grid(true)\n", "xlabel(\"integration step\")\n", "ylabel(\"absolute error\")\n", "title(\"Error of trapezoidal rule\")\n", "legend();" ] }, { "cell_type": "markdown", "id": "f388c852-0232-42dc-b34b-06a5870bfea4", "metadata": {}, "source": [ "## Richardson extrapolation" ] }, { "cell_type": "markdown", "id": "4b074d95-a332-468a-a91e-1ca3420792cb", "metadata": {}, "source": [ "Recall that \n", "$$\\int_a^b f(x) \\, \\mathrm{d}x = T(h, n) + \\alpha h^2 + O(h^p),$$ where $T(h, n)$ is given by Eq. (2).\n", "\n", "Next, we write down the trapezoidal rule with step $2h$. It will involve\n", "only the nodes at $x_1$, $x_3$, $x_5$, $\\ldots$, $x_{n-2}$, $x_n$:\n", "\n", "$$\\int_a^b f(x) \\, \\mathrm{d}x = h [ f(x_1) + 2f(x_3) + 2f(x_5) + \\ldots + 2f(x_{n-2}) + f(x_n) ] + \\alpha (2h)^2 + O(h^p) \n", "= T(2h) + 4\\alpha h^2 + O(h^p) . \\tag{3}$$\n", "\n", "The leading error in Eq. (3) is 4 times the error in Eq. (2), since we doubled the step size. So we multiply both sides of Eq. (2) by 4,\n", "\n", "$$\n", "4\\int_a^b f(x) \\, \\mathrm{d}x = 4T(h) + 4\\alpha h^2 + O(h^p), \\tag{4}\n", "$$\n", "and subtract Eq. (3) from Eq. (4), intending to\n", "cancel out the leading error term. The result that we get is the\n", "*Simpson's formula*:\n", "\n", "$$\n", "\\int_a^b f(x) \\, \\mathrm{d}x = \\frac{1}{3}\\left(4T(h) - T(2h) \\right) + O(h^p) .\n", "$$" ] }, { "cell_type": "markdown", "id": "a64dac0c-4520-4427-b9b5-6ab62e41d452", "metadata": {}, "source": [ "### Simpson's formula" ] }, { "cell_type": "code", "execution_count": null, "id": "9e499f89-d6af-4519-bf32-79242640a90b", "metadata": {}, "outputs": [], "source": [ "\"\"\"\n", " ans = mytrapezoids2(f, h)\n", "\n", "Numerically evaluate the integral int_a^b fun(x) dx using the\n", "trapezoidal rule: I = h/2*(f[1] + 2f[2] + ... + 2f[n-1] + f[n]),\n", "where h = (b - a)/(n - 1), x_i = a + (i - 1)*h, f[i] = fun(x_i).\n", "\"\"\"\n", "function mytrapezoids2(fvec, h)\n", " s1 = fvec[1] + fvec[end]\n", " s2 = sum(fvec)\n", " return (2*s2 - s1)*h/2\n", "end" ] }, { "cell_type": "code", "execution_count": null, "id": "4c1f4b14-f1c6-4d00-8077-7fae9793eb96", "metadata": {}, "outputs": [], "source": [ "\"\"\"\n", " ans = mysimpsons2(fvec, h)\n", "\n", "Numerically evaluate the integral int_a^b fun(x) dx using the\n", "Simpson's rule, where h = (b - a)/(n - 1), x_i = a + (i - 1)*h, fvec[i] = fun(x_i).\n", "\"\"\"\n", "function mysimpsons2(fvec, h)\n", " n = length(fvec)\n", " t1 = mytrapezoids2(fvec, h)\n", " t2 = mytrapezoids2(fvec[1:2:n], 2*h)\n", " return (4*t1 - t2)/3\n", "end" ] }, { "cell_type": "code", "execution_count": null, "id": "5b7f54dd-c5b6-4788-8d89-a8a53167405a", "metadata": {}, "outputs": [], "source": [ "# Simple test 1\n", "fun1(x) = x\n", "a = -1.0\n", "b = 1.0\n", "np = 21\n", "exact = 0.0\n", "xx = range(a, b, np)\n", "h = xx[2] - xx[1]\n", "fvec = fun1.(xx)\n", "abs(mysimpsons2(fvec, h) - exact) " ] }, { "cell_type": "code", "execution_count": null, "id": "2447524e-fe7c-4ce2-b85a-a7d020b695a5", "metadata": {}, "outputs": [], "source": [ "# Simple test 2\n", "fun2(x) = exp(x)\n", "a = 0.0\n", "b = 1.0\n", "np = 201\n", "exact = exp(1.0) - 1.0\n", "xx = range(a, b, np)\n", "h = xx[2] - xx[1]\n", "fvec = fun2.(xx)\n", "abs(mysimpsons2(fvec, h) - exact) " ] }, { "cell_type": "markdown", "id": "85da9da2-cbfc-4ac9-94e4-3c7b9d813166", "metadata": {}, "source": [ "### Error of the Simpson's formula" ] }, { "cell_type": "markdown", "id": "8812a4df-d5f6-4d96-9d48-c02f78a6af8e", "metadata": {}, "source": [ "An analysis of the error of Simpson's rule vs integration step h" ] }, { "cell_type": "code", "execution_count": null, "id": "497e904c-fa1e-4ce3-a307-9df2e1502cde", "metadata": {}, "outputs": [], "source": [ "fun0(x) = sin(x)\n", "a = 0.0\n", "b = pi\n", "exact = 2.0\n", "ndp = 10\n", "errs = zeros(ndp)\n", "hh = zeros(ndp);" ] }, { "cell_type": "code", "execution_count": null, "id": "fa9f5752-0a00-4f5c-8e8e-c4c3d4307922", "metadata": {}, "outputs": [], "source": [ "for k = 1:ndp \n", " np = 2^k + 1\n", " xx = range(a, b, np)\n", " h = xx[2] - xx[1]\n", " hh[k] = h\n", " fun0vec = fun0.(xx)\n", " errs[k] = abs(mysimpsons2(fun0vec, h) - exact)\n", "end" ] }, { "cell_type": "code", "execution_count": null, "id": "a76abeeb-c237-4233-adf0-721eee76a1bd", "metadata": {}, "outputs": [], "source": [ "loglog(hh, errs, label=\"numerical experiment\", marker=\".\")\n", "loglog(hh, hh .^ 4, linestyle=\"dashed\", label=L\"hh^4\")\n", "grid(true)\n", "xlabel(\"integration step\")\n", "ylabel(\"absolute error\")\n", "title(\"Error of Simpson's rule\")\n", "legend();" ] }, { "cell_type": "code", "execution_count": null, "id": "40b1ec1a-bc2e-4c19-be66-ead5a219825e", "metadata": {}, "outputs": [], "source": [ "loglog(hh, errs, label=\"Simpson's\", marker=\".\")\n", "loglog(hh, errt, label=\"Trapezoidal\", marker=\".\")\n", "loglog(hh, hh .^ 4 / 1000, linestyle=\"dashed\", label=L\"hh^4\")\n", "loglog(hh, hh .^ 2, linestyle=\"dotted\", label=L\"hh^2\")\n", "grid(true)\n", "xlabel(\"integration step\")\n", "ylabel(\"absolute error\")\n", "title(\"Errors of Simpson's and trapezoidal rules\")\n", "legend();" ] }, { "cell_type": "code", "execution_count": null, "id": "ff45c07b-879c-486e-bb84-9f8a2033e515", "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 }