{
"cells": [
{
"cell_type": "markdown",
"id": "c742c64c",
"metadata": {},
"source": [
"# Quadrature Rules"
]
},
{
"cell_type": "markdown",
"id": "4c1f2ba2",
"metadata": {},
"source": [
"This notebook introduces rules to approximate definite integrals."
]
},
{
"cell_type": "markdown",
"id": "8393d2dc",
"metadata": {},
"source": [
"## The plot of an example function"
]
},
{
"cell_type": "code",
"execution_count": 1,
"id": "36efa60e",
"metadata": {},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"\n",
"(process:82088): GLib-GIO-WARNING **: 16:50:59.652: Unexpectedly, UWP app `Clipchamp.Clipchamp_2.5.1.0_neutral__yxz26nhyzhsrt' (AUMId `Clipchamp.Clipchamp_yxz26nhyzhsrt!App') supports 41 extensions but has no verbs\n",
"\n",
"(process:82088): GLib-GIO-WARNING **: 16:50:59.652: Unexpectedly, UWP app `Evernote.Evernote_10.46.9.0_x64__q4d96b2w5wcc2' (AUMId `Evernote.Evernote_q4d96b2w5wcc2!Evernote') supports 1 extensions but has no verbs\n"
]
}
],
"source": [
"using Plots"
]
},
{
"cell_type": "code",
"execution_count": 2,
"id": "662de2fc",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"f (generic function with 1 method)"
]
},
"execution_count": 2,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"f(x) = (x-1)^3 + 0.1*x^2 - x + 2"
]
},
{
"cell_type": "code",
"execution_count": 3,
"id": "d9279ab4",
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 3,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"r = -0.35:0.01:2.2\n",
"y = [f(t) for t in r]\n",
"plot(r,y,label=\"y=f(x)\")"
]
},
{
"cell_type": "markdown",
"id": "89193a20",
"metadata": {},
"source": [
"## The Midpoint Rule"
]
},
{
"cell_type": "code",
"execution_count": 4,
"id": "5add65f7",
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 4,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"plot(r,y,label=\"y=f(x)\")\n",
"intrange = 0:0.01:2\n",
"mf0 = [0 for t in intrange]\n",
"mf1 = [f(1) for t in intrange]\n",
"plot!(intrange, mf0, line=(:red), label=\"\")\n",
"plot!(intrange, mf1, line=(:red), label=\"\")\n",
"plot!([0.0; 0.0], [0.0; f(1)], line=(:red), label=\"\")\n",
"plot!([2.0; 2.0], [0.0; f(1)], line=(:red), label=\"\")\n",
"plot!([1.0; 1.0], [0.0; f(1)], line=(:black, :dot), label=\"\")\n",
"scatter!([1.0; 1.0],[f(1); f(1)],marker=(:blue),label=\"f(1)\")\n",
"# savefig(\"figmidpointrule.png\")"
]
},
{
"cell_type": "markdown",
"id": "aac88b3b",
"metadata": {},
"source": [
"## The Composite Midpoint Rule"
]
},
{
"cell_type": "code",
"execution_count": 5,
"id": "fd091c64",
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 5,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"plot(r,y,label=\"y=f(x)\")\n",
"intrange = 0:0.01:2\n",
"mf0 = [0 for t in intrange]\n",
"intrange1 = 0.0:0.01:0.5\n",
"intrange2 = 0.5:0.01:1.0\n",
"intrange3 = 1.0:0.01:1.5\n",
"intrange4 = 1.5:0.01:2.0\n",
"mf1 = [f(0.25) for t in intrange1]\n",
"mf2 = [f(0.75) for t in intrange2]\n",
"mf3 = [f(1.25) for t in intrange3]\n",
"mf4 = [f(1.75) for t in intrange4]\n",
"plot!(intrange, mf0, line=(:red), label=\"\")\n",
"plot!(intrange1, mf1, line=(:red), label=\"\")\n",
"plot!(intrange2, mf2, line=(:red), label=\"\")\n",
"plot!(intrange3, mf3, line=(:red), label=\"\")\n",
"plot!(intrange4, mf4, line=(:red), label=\"\")\n",
"plot!([0.0; 0.0], [0.0; f(0.25)], line=(:red), label=\"\")\n",
"plot!([0.5; 0.5], [0.0; f(0.25)], line=(:red), label=\"\")\n",
"plot!([1.0; 1.0], [0.0; f(0.75)], line=(:red), label=\"\")\n",
"plot!([1.5; 1.5], [0.0; f(1.75)], line=(:red), label=\"\")\n",
"plot!([2.0; 2.0], [0.0; f(1.75)], line=(:red), label=\"\")\n",
"plot!([0.25; 0.25], [0.0; f(0.25)], line=(:black, :dot), label=\"\")\n",
"plot!([0.75; 0.75], [0.0; f(0.75)], line=(:black, :dot), label=\"\")\n",
"plot!([1.25; 1.25], [0.0; f(1.25)], line=(:black, :dot), label=\"\")\n",
"plot!([1.75; 1.75], [0.0; f(1.75)], line=(:black, :dot), label=\"\")\n",
"scatter!([0.25; 0.25],[f(0.25); f(0.25)],marker=(:blue),label=\"f(0.25)\")\n",
"scatter!([0.75; 0.75],[f(0.75); f(0.75)],marker=(:green),label=\"f(0.75)\")\n",
"scatter!([1.25; 1.25],[f(1.25); f(1.25)],marker=(:yellow),label=\"f(1.25)\")\n",
"scatter!([1.75; 1.75],[f(1.75); f(1.75)],marker=(:brown),label=\"f(1.75)\")\n",
"# savefig(\"figmidpointcomposite.png\")"
]
},
{
"cell_type": "markdown",
"id": "cba51de8",
"metadata": {},
"source": [
"## The Trapezoidal Rule"
]
},
{
"cell_type": "code",
"execution_count": 6,
"id": "bf6b00a3",
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 6,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"plot(r,y,label=\"y=f(x)\")\n",
"intrange = 0:0.01:2\n",
"mf0 = [0 for t in intrange]\n",
"plot!(intrange, mf0, line=(:red), label=\"\")\n",
"yT(x) = f(0) + (f(2)-f(0))*x/2\n",
"my = [yT(t) for t in intrange]\n",
"plot!(intrange, my, line=(:red), label=\"\")\n",
"r0 = [0.0; 0.0]\n",
"f0 = [0.0; f(1)]\n",
"plot!([0.0; 0.0], [0.0; f(0)], line=(:red), label=\"\")\n",
"plot!([2.0; 2.0], [0.0; f(2)], line=(:red), label=\"\")\n",
"plot!([0.0; 0.0], [0.0; f(0)], line=(:black, :dot), label=\"\")\n",
"plot!([2.0; 2.0], [0.0; f(2)], line=(:black, :dot), label=\"\")\n",
"scatter!([0.0; 0.0],[f(0); f(0)],marker=(:blue),label=\"f(0)\")\n",
"scatter!([2.0; 2.0],[f(2); f(2)],marker=(:green),label=\"f(2)\")\n",
"# savefig(\"figtrapezoidalrule.png\")"
]
},
{
"cell_type": "markdown",
"id": "88f1425e",
"metadata": {},
"source": [
"## Simpson's Rule"
]
},
{
"cell_type": "markdown",
"id": "32b6dae2",
"metadata": {},
"source": [
"Simpson's rule to approximate $\\int_a^b f(x) dx$ evaluates $f$ at $a$, $(a+b)/2$ and $b$. We will use SymPy to determine the weights $w_0$, $w_1$, and $w_2$ so the rule\n",
"\n",
"$$\n",
" w_0 f(a) + w_1 f\\left( \\frac{a+b}{2} \\right) + w_2 f(b) \\approx \\int_a^b f(x) dx\n",
"$$\n",
"has a degree of precision as high as possible. This means that the rule is exact for the first three basis functions: 1, $x$, and $x^2$."
]
},
{
"cell_type": "code",
"execution_count": 7,
"id": "ab75b565",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"equation 1 : w0 + w1 + w2\n",
"equation 2 : a*w0 + b*w2 + w1*(a + b)/2\n",
"equation 3 : a^2*w0 + b^2*w2 + w1*(a + b)^2/4\n"
]
}
],
"source": [
"using SymPy\n",
"a,b,w0,w1,w2 = Sym(\"a,b,w0,w1,w2\")\n",
"eq1 = w0 + w1 + w2\n",
"eq2 = a*w0 + (a+b)*w1/2 + b*w2\n",
"eq3 = a^2*w0 + (a+b)^2*w1/4 + b^2*w2\n",
"println(\"equation 1 : \", eq1)\n",
"println(\"equation 2 : \", eq2)\n",
"println(\"equation 3 : \", eq3)"
]
},
{
"cell_type": "code",
"execution_count": 8,
"id": "58b8a162",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"Dict{Any, Any} with 3 entries:\n",
" w0 => -a/6 + b/6\n",
" w1 => -2*a/3 + 2*b/3\n",
" w2 => -a/6 + b/6"
]
},
"execution_count": 8,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"solve([eq1 - b+a,\n",
" eq2 - b^2/2 + a^2/2,\n",
" eq3 - b^3/3 + a^3/3], [w0,w1,w2])"
]
},
{
"cell_type": "markdown",
"id": "e010e9e2",
"metadata": {},
"source": [
"Then Simpson's rule is\n",
"\n",
"$$\n",
" \\left( \\frac{b-a}{6} \\right) f(a) + 2 \\left( \\frac{b-a}{3} \\right) f\\left( \\frac{a+b}{2} \\right)\n",
" + \\left( \\frac{b-a}{6} \\right) f(b) \\approx \\int_a^b f(x) dx.\n",
"$$"
]
},
{
"cell_type": "markdown",
"id": "a2956ff0",
"metadata": {},
"source": [
"## The Composite Trapezoidal Rule"
]
},
{
"cell_type": "code",
"execution_count": 9,
"id": "cd00b89f",
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 9,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"plot(r,y,label=\"y=f(x)\")\n",
"intrange = 0:0.01:2\n",
"mf0 = [0 for t in intrange]\n",
"plot!(intrange, mf0, line=(:red), label=\"\")\n",
"yT1(x) = f(0) + (f(0.5)-f(0))*x/0.5\n",
"intrange1 = 0:0.01:0.5\n",
"my1 = [yT1(t) for t in intrange1]\n",
"plot!(intrange1, my1, line=(:red), label=\"\")\n",
"yT2(x) = f(0.5) + (f(1.0)-f(0.5))*(x-0.5)/0.5\n",
"intrange2 = 0.5:0.01:1.0\n",
"my2 = [yT2(t) for t in intrange2]\n",
"plot!(intrange2, my2, line=(:red), label=\"\")\n",
"yT3(x) = f(1.0) + (f(1.5)-f(1.0))*(x-1.0)/0.5\n",
"intrange3 = 1.0:0.01:1.5\n",
"my3 = [yT3(t) for t in intrange3]\n",
"plot!(intrange3, my3, line=(:red), label=\"\")\n",
"yT4(x) = f(1.5) + (f(2.0)-f(1.5))*(x-1.5)/0.5\n",
"intrange4 = 1.5:0.01:2.0\n",
"my4 = [yT4(t) for t in intrange4]\n",
"plot!(intrange4, my4, line=(:red), label=\"\")\n",
"plot!([0.0; 0.0], [0.0; f(0)], line=(:red), label=\"\")\n",
"plot!([0.5; 0.5], [0.0; f(0.5)], line=(:red), label=\"\")\n",
"plot!([1.0; 1.0], [0.0; f(1.0)], line=(:red), label=\"\")\n",
"plot!([1.5; 1.5], [0.0; f(1.5)], line=(:red), label=\"\")\n",
"plot!([2.0; 2.0], [0.0; f(2.0)], line=(:red), label=\"\")\n",
"plot!([0.0; 0.0], [0.0; f(0)], line=(:black, :dot), label=\"\")\n",
"plot!([0.5; 0.5], [0.0; f(0.5)], line=(:black, :dot), label=\"\")\n",
"plot!([1.0; 1.0], [0.0; f(1.0)], line=(:black, :dot), label=\"\")\n",
"plot!([1.5; 1.5], [0.0; f(1.5)], line=(:black, :dot), label=\"\")\n",
"plot!([2.0; 2.0], [0.0; f(2.0)], line=(:black, :dot), label=\"\")\n",
"scatter!([0.0; 0.0],[f(0); f(0)],marker=(:blue),label=\"f(0)\")\n",
"scatter!([0.5; 0.5],[f(0.5); f(0.5)],marker=(:green),label=\"f(0.5)\")\n",
"scatter!([1.0; 1.0],[f(1.0); f(1.0)],marker=(:yellow),label=\"f(1.0)\")\n",
"scatter!([1.5; 1.5],[f(1.5); f(1.5)],marker=(:brown),label=\"f(1.5)\")\n",
"scatter!([2.0; 2.0],[f(2.0); f(2.0)],marker=(:purple),label=\"f(2.0)\")\n",
"# savefig(\"figtrapezoidalcomposite.png\")"
]
},
{
"cell_type": "markdown",
"id": "0598a09f",
"metadata": {},
"source": [
"The Julia function below applies the composite Trapezoidal rule to some function over an interval."
]
},
{
"cell_type": "code",
"execution_count": 10,
"id": "b01cb592",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"comptrap"
]
},
"execution_count": 10,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"\"\"\"\n",
"Applies the composite trapezoidal rule to approximate the integral\n",
"of the function f over the interval [a,b], using n subintervals.\n",
"\n",
"Example: t = comptrap(cos,0.0,pi/2,100)\n",
"\"\"\"\n",
"function comptrap(f::Function,a::Float64,b::Float64,n::Int64)\n",
" h = (b-a)/n\n",
" t = (f(a) + f(b))/2\n",
" for i = 1:n-1\n",
" t = t + f(a+i*h)\n",
" end\n",
" t = t*h\n",
" return t\n",
"end"
]
},
{
"cell_type": "markdown",
"id": "6d939329",
"metadata": {},
"source": [
"Let us run the example."
]
},
{
"cell_type": "code",
"execution_count": 11,
"id": "4065cb67",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"0.9999794382396076"
]
},
"execution_count": 11,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"t = comptrap(cos,0.0,pi/2,100)"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Julia 1.8.0",
"language": "julia",
"name": "julia-1.8"
},
"language_info": {
"file_extension": ".jl",
"mimetype": "application/julia",
"name": "julia",
"version": "1.8.0"
},
"latex_envs": {
"LaTeX_envs_menu_present": true,
"autocomplete": true,
"bibliofile": "biblio.bib",
"cite_by": "apalike",
"current_citInitial": 1,
"eqLabelWithNumbers": true,
"eqNumInitial": 1,
"hotkeys": {
"equation": "Ctrl-E",
"itemize": "Ctrl-I"
},
"labels_anchors": false,
"latex_user_defs": false,
"report_style_numbering": false,
"user_envs_cfg": false
}
},
"nbformat": 4,
"nbformat_minor": 5
}