\n",
"Using calculus we see that \n",
"$$\\int_{-\\frac{\\pi}{2}}^\\frac{\\pi}{2} \\cos x~dx=\n",
"[\\sin x]_{-\\frac{\\pi}{2}}^\\frac{\\pi}{2} =\n",
"\\sin(\\frac{\\pi}{2})-\\sin(\\frac{-\\pi}{2})=1-(-1)=2.$$\n",
"The *absolute error* is the absolute value of the difference between an approximation and the actual value, so:\n",
"

"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"0.09439510239319526"
]
},
"execution_count": 7,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"absolute_error = abs(approx - 2)\n",
"absolute_error"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Degree of precision\n",
"\n",
"An integration rule over an interval $[a,b]$ has *degree of precision* $m$ if it is exact (zero error) for all polynomials of degree at most $m$, but is not exact for $x^{m+1}$. "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"*Remark:* To check exactness of all polynomials of degree up to $m$, it is enough to check the polynomials\n",
"$$1, x, x^2, \\ldots, x^m.$$\n",
"This is because of our choice of approximation as a weighted sum, $$A(f) = \\sum_{k=0}^N c_k f(x_k).$$ \n",
"This expression is linear in $f$, meaning that:\n",
"* If $f$ is a function and $s \\in {\\mathbb R}$ is a scalar, then \n",
"$A(s f)=s A(f)$ where $sf$ denotes the function sending $x$ to $s f(x)$.\n",
"* If $f$ and $g$ are functions, then $A(f+g)=A(f)+A(g)$, where $f+g$ denotes the function sending $x$ to $f(x)+g(x)$.\n",
"The integral is also linear in the sense above. Using the two rules above, you can show (for example) that if\n",
"$$A(1) = \\int_a^b 1~dx \n",
"\\quad \\text{and} \\quad \n",
"A(x) = \\int_a^b x~dx,$$\n",
"then $A$ is exact for any linear function, i.e., $A(mx+c)=\\int_a^b mx+c~dx$ for all $m,c \\in {\\mathbb R}.$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
" **Result.** The Midpoint Rule and Trapezoid Rule have degree of precision $1$.

\n",
" \n",
"**Proof for the midpoint rule:**\n",
"Let $A$ denote the midpoint rule applied to the interval $[a,b]$. \n",
"We need to show that $A$ is exact for the function $1$, and the function $x$, but not for $x^2$.\n",
"\n",
"First we check it for the constant function $1$. Observe $\\int_a^b 1~dx=b-a$ and \n",
"$$A(1)=(b-a)\\cdot 1=b-a.$$\n",
"So $A$ is exact for $1$.\n",
"\n",
"Now we check it for the function $x$. Observe that \n",
"$$\\int_a^b x~dx=\\left[\\frac{1}{2}x^2\\right]_a^b = \\frac{1}{2}(b^2-a^2).$$\n",
"We also have\n",
"$$A(x) = (b-a) \\cdot \\frac{b-a}{2}=\\frac{1}{2}(b^2-a^2).$$\n",
"Thus, $A$ is exact for $x$.\n",
"\n",
"To see that $A$ is not exact for $x^2$, observe that \n",
"$$\\int_a^b x^2~dx=\\left[\\frac{1}{3}x^3\\right]_a^b = \\frac{1}{3}(b^3-a^3)$$\n",
"and\n",
"$$A(x^2) = (b-a) \\cdot \\left(\\frac{b-a}{2}\\right)^2=\\frac{1}{4}(b-a)^3.$$\n",
"**QED**"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We will not prove the trapezoid case. It should be clear geometrically."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
" **Result.** The Simpson's Rule has degree of precision $3$.

\n",
"\n",
"This is explained in TAK on page 54, so we will not discuss it here."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Problem\n",
"\n",
"Derive a rule which has degree of precision $1$ using the points $1$ and $3$ in the interval $[0,5]$."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"**Solution.**\n",
"\n",
"A rule using points $1$ and $3$ must have the form \n",
"$$A(f) = c_1 f(1) + c_3 f(3).$$\n",
"Since we want our rule to have degree $1$, it must satisfy \n",
"$$A(1)=\\int_0^5 1~dx =5 \n",
"\\quad \\text{and} \\quad \n",
"A(x)=\\int_0^5 x~dx =\\frac{25}{2}.$$\n",
"Observe that by definition, we have\n",
"$$A(1)=c_1 + c_3 \n",
"\\quad \\text{and} \\quad \n",
"A(x) = c_1 + 3 c_3.$$\n",
"So, we must solve the equations\n",
"$$c_1+c_3 = 5\n",
"\\quad \\text{and} \\quad \n",
"c_1 + 3 c_3=\\frac{25}{2}.$$\n",
"Subtracting the first equation from the second yields\n",
"$$2 c_3 = \\frac{15}{2} \\quad \\text{so} \\quad c_3=\\frac{15}{4}.$$\n",
"Then plugging this value of $c_3$ into one of the equations and solving yields\n",
"$$c_1=\\frac{5}{4}.$$\n",
"Thus, our approximation rule has the form \n",
"$$A(f) = \\frac{5}{4} f(1) + \\frac{15}{4} f(3).$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Below we check that our formula works as expected:"
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {},
"outputs": [],
"source": [
"def A(f):\n",
" return 5/4*f(1) + 15/4*f(3)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Here we expect to get $\\int_0^5 1~dx=5$:"
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"5.0"
]
},
"execution_count": 9,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"A(lambda x: 1)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Here we expect to get $\\int_0^5 x~dx=12\\frac{1}{2}$:"
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"12.5"
]
},
"execution_count": 10,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"A(lambda x: x)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We can see that our formula is not exact in degree $2$ since\n",
"$\\int_0^5 x^2~dx=\\frac{125}{3}=41\\frac{2}{3}$, but:"
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"35.0"
]
},
"execution_count": 11,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"A(lambda x: x**2)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Gaussian Quadrature\n",
"\n",
"The idea of Gaussian Quadrature is to find the approximation with the highest degree of precision among all approximations using $N$ points. It is know that you can always find an approximation with degree $2N-1$.\n",
"\n",
"The book works out the formula for the Gaussian Quadrature using $2$ and $3$ points in the interval $[-1,1]$. The two point formula is\n",
"$$G_2(f) = f\\left(\\frac{1}{\\sqrt{3}}\\right) + f\\left(\\frac{-1}{\\sqrt{3}}\\right),$$\n",
"which has degree of precision $3$.\n",
"The three-point formula is\n",
"$$G_3(f) = \\frac{5}{9}\\,f\\left(\\sqrt{\\frac{3}{5}}\\right) + \\frac{8}{9}\\,f(0) + \\frac{5}{9}\\,f\\left(-\\sqrt{\\frac{3}{5}}\\right),$$\n",
"which has degree of precision $5$.\n",
"These are derived in example 9 (starting on page 55 of TAK)."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Changes of coordinates\n",
"\n",
"\n", "Recall the change of coordinates formula for integration. If $f$ is integrable and $g$ is differentiable, then \n", "$$\\int_a^b f\\big(g(x)\\big)g'(x)~dx = \\int_{g(a)}^{g(b)} f(x)~dx.$$\n", "We will apply this to the $3$-point Gaussian Quadrature to derive a formula which holds over any interval.\n", "

\n", "\n", "\n", "We can apply this to move an integral defined over say $[-1,1]$ to any other interval. If we want to consider an integral defined over $[c,d]$ we can take \n", "$$g(x)=\\frac{d-c}{2} x + \\frac{c+d}{2}.$$\n", "This is a linear function, thus if $f(x)$ is a polynomial $f\\big(g(x)\\big)g'(x)$ is a polynomial of the same degree. Applying our change of coordinates formula, and pulling out the constant $g'(x)=\\frac{d-c}{2}$ yields\n", "$$\\frac{d-c}{2} \\int_{-1}^1 f\\big(g(x)\\big)~dx = \\int_{c}^{d} f(x)~dx.$$\n", "The integral on the left is over the interval $[-1,1]$ so our $G_3$ will be exact provided $f\\big(g(x)\\big)$ is a polynomial of degree 5 or less. Thus for polynomials $f$ of degree $5$ or less, we have\n", "$$\\frac{d-c}{2} \\left[\\frac{5}{9}\\,f \\circ g\\left(\\sqrt{\\frac{3}{5}}\\right) + \\frac{8}{9}\\,f \\circ g(0) + \\frac{5}{9}\\,f\\circ g\\left(-\\sqrt{\\frac{3}{5}}\\right)\\right] = \\int_{c}^{d} f(x)~dx.$$\n", "We find it easiest to leave our expression like this, but we could expand $g(\\star)$ to get points, so this is still in the form of an integration rule. Also it still has degree of precision $5$. This is a special case of a more general rule:\n", "

\n", "\n", "\n",
"**Change of coordinates formula for quadrature rules.**\n",
"Suppose $A(f)=\\sum_{k=0}^N c_k f(x_k)$ is an quadrature rule with degree of precision $m$ over $[-1,1]$.\n",
"Let $[c,d]$ be another interval and define $g(x)=\\frac{d-c}{2} x + \\frac{c+d}{2}.$\n",
"Let $y_k=g(x_k)$ and $d_k = \\frac{d-c}{2} c_k$. Then the quadrature rule\n",
"$B(f)=\\sum_{k=0}^N d_k f(y_k)$ has degree of precision $m$ over $[c,d]$.\n",
"

\n",
"\n",
"\n", "We apply this below to the interval from $[0,\\pi]$ to estimate an integral.\n", "

" ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [], "source": [ "c = 0\n", "d = m.pi\n", "\n", "# Linear change of coordinates [-1,1] -> [c,d]\n", "g = lambda x: (d-c)/2*x + (c+d)/2\n", "\n", "# Points used in Gaussian Quadrature over [-1,1]\n", "x1 = -m.sqrt(3)/m.sqrt(5)\n", "x2 = 0\n", "x3 = m.sqrt(3)/m.sqrt(5)\n", "\n", "# Weights used in Gaussian Quadrature over [-1,1]\n", "c1 = 5/9\n", "c2 = 8/9\n", "c3 = 5/9\n", "\n", "# Points used in Gaussian Quadrature over [c,d]\n", "y1 = g(x1)\n", "y2 = g(x2)\n", "y3 = g(x3)\n", "\n", "# Weights used in Gaussian Quadrature over [c,d]\n", "d1 = (d-c)/2 * c1\n", "d2 = (d-c)/2 * c2\n", "d3 = (d-c)/2 * c3\n", "\n", "# Gaussian Quadrature formula over [c,d]\n", "A = lambda f: d1*f(y1) + d2*f(y2) + d3*f(y3)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Our $A$ should be have degree of precision $5$. We will check that. \n", "\n", "Below we compare $A(x^k)$ and $\\int_0^\\pi x^k~dx=\\frac{1}{k+1}\\pi^{k+1}$ for values of $k \\leq 6$." ] }, { "cell_type": "code", "execution_count": 13, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "For x^0, our approximation is 3.141592653589793 and the exact value is 3.141592653589793.\n", "The error was 0.0\n", "For x^1, our approximation is 4.934802200544679 and the exact value is 4.934802200544679.\n", "The error was 0.0\n", "For x^2, our approximation is 10.335425560099939 and the exact value is 10.335425560099939.\n", "The error was 0.0\n", "For x^3, our approximation is 24.352272758500604 and the exact value is 24.352272758500604.\n", "The error was 0.0\n", "For x^4, our approximation is 61.20393695705627 and the exact value is 61.203936957056285.\n", "The error was -1.4210854715202004e-14\n", "For x^5, our approximation is 160.23153226255067 and the exact value is 160.2315322625507.\n", "The error was -2.842170943040401e-14\n", "For x^6, our approximation is 430.39178495819266 and the exact value is 431.4704611109702.\n", "The error was -1.0786761527775184\n" ] } ], "source": [ "for k in range(7):\n", " approx = A(lambda x: x**k)\n", " exact = m.pi**(k+1) / (k+1)\n", " print(\"For x^{}, our approximation is {} and the exact value is {}.\".format(k, approx, exact))\n", " print(\"The error was {}\".format(approx-exact))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "It seems clear that it was exact up to degree $5$ (with non-zero values of the order $10^{-14}$ probably coming from errors in floating point arithmetic.) The degree $6$ error of about $1$ indicates that our formula is not exact for polynomials of degree $6$, which we expect." ] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "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.6.9" } }, "nbformat": 4, "nbformat_minor": 2 }