{
"cells": [
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
" \n",
"\n"
],
"text/plain": [
""
]
},
"execution_count": 3,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"from IPython.core.display import HTML\n",
"def css_styling():\n",
" try:\n",
" with open(\"tma4125.css\", \"r\") as f:\n",
" styles = f.read()\n",
" return HTML(styles)\n",
" except FileNotFoundError:\n",
" pass #Do nothing\n",
"\n",
"# Comment out next line and execute this cell to restore the default notebook style \n",
"css_styling()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Exercises 3\n",
"\n",
"Submission deadline: **Feb 15 2021 at 12:00 (noon)**\n",
"\n",
"In this exercise set we will construct and analyze quadrature rules. For guidance on quadrature rules, please read the lecture notes. Make sure to run the code below to get all the important modules, and to make the plots look nice."
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [],
"source": [
"%matplotlib inline\n",
"\n",
"from numpy import *\n",
"from matplotlib.pyplot import *\n",
"from math import factorial\n",
"newparams = {'figure.figsize': (8.0, 4.0), 'axes.grid': True,\n",
" 'lines.markersize': 8, 'lines.linewidth': 2,\n",
" 'font.size': 14}\n",
"rcParams.update(newparams)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## 1) Analyzing the composite Simpson's rule\n",
"\n",
"Simpson's rule is defined as\n",
"$$\n",
"S[f](x_{i-1}, x_i) = \\frac{h}{6}(f(x_{i-1} + 4f(x_{i-1/2}) + f(x_i))\n",
"$$\n",
"where $h = x_i - x_{i-1}$ and $x_{i-1/2} = \\frac{x_{i-1}+x_i}{2}.$\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"**a)**\n",
"\n",
"Show that the resulting composite Simpson's rule is given by\n",
"\n",
"\\begin{align*}\n",
"\\int_a^b f {\\,\\mathrm{d}x} \\approx \\mathrm{CSR}[f]({[x_{i-1}, x_i]}_{i=1}^{m})\n",
"&= \n",
"\\tfrac{h}{6}\n",
"[\n",
"f(x_0)\n",
"+ 4f(x_{x_{1/2}}) + 2f(x_1) \n",
"+ 4f(x_{3/2}) + 2f(x_2)\n",
"+ \\ldots\n",
"\\\\ \n",
"&\\qquad+2 f(x_{m-1})\n",
"+\n",
"4f(x_{x_{m-1/2}}) \n",
"+f(x_m)\n",
"].\n",
"\\end{align*}"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"**b)**\n",
"Implement the composite Simpson's rule.\n",
"Use this function to compute an approximate value of integral\n",
"\n",
"$$\n",
"I(0,1) = \\int_0^1 \\cos\\left(\\frac{\\pi}{2}x\\right) = \\frac{2}{\\pi} = 0.636619\\dotsc.\n",
"$$\n",
"\n",
"for $m = 4, 8, 16, 32, 64$ corresponding to\n",
"$ h = 2^{-2}, 2^{-3}, 2^{-4}, 2^{-5}, 2^{-6}$.\n",
"Tabulate the corresponding quadrature errors $|I(0,1) - Q(0,1)|$. Plot the errors against $h$.\n",
"What do you observe? How does it compare to the composite trapezoidal rule?"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Doubling the number of sub-intervals decreases the error by a factor of 16. This suggests that the error of the composite Simpson rule is $Cm^{-4} = Ch^4$."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"**c)**\n",
"Recall that the error of Simpson's rule on a single interval is given by\n",
"\n",
"$$\n",
"|I[f](a,b) - S[f](a,b)| = - \\frac{(b-a)^5}{2880} f^{(4)}(\\xi)\n",
"$$\n",
"for some $\\xi \\in [a,b]$.\n",
"\n",
"Use this to show that the error of the composite Simpson rule can be bounded by\n",
"\n",
"\\begin{equation}\n",
"|I[f]-\\mathrm{CSR}[f]|\n",
"\\leqslant\n",
"\\dfrac{M_4}{2880} \\dfrac{(b-a)^5}{m^4}\n",
"=\n",
"\\dfrac{M_4}{2880}\n",
"h^4(b-a)\n",
"\\label{quad:thm:csr_est} \\tag{3}\n",
"\\end{equation}\n",
"\n",
"where $M_4 = \\max_{\\xi\\in[a,b]} |f^{(4)}(\\xi)|$."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## 2) Gaussian Quadrature\n",
"\n",
"In this exercise we will construct a Gaussian quadrature rule with 3 nodes. We will take it step by step, so don't worry if you do not feel like an expert on Gaussian quadrature. \n",
"\n",
"To make your life easy, we will use the [sympy](https://docs.sympy.org/latest/index.html#) python module for symbolic mathematics to perform tasks such as (symbolic) integration and root finding of low order polynoms. \n",
"In particular look at [integrate](https://docs.sympy.org/latest/modules/integrals/integrals.html) and [solve](https://docs.sympy.org/latest/modules/solvers/solvers.html) submodules."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"\n",
"The first step in constructing a Gaussian quadrature is finding the correct orthogonal polynomial. The nodes of the quadrature rule will be the roots of some polynomial. Since we are looking for 3 nodes, this means that the polynomial should have 3 roots, and hence we are looking for a third-order polynomial.\n",
" \n",
"The polynomial, call it $p_3$, should be orthogonal on the interval $[0,3]$ to all polynomials of order 2 or less. We now create this polynomial.\n",
"\n",
"Start with the 4 polynomial \"basis\" functions \n",
"$$\n",
"\\phi_0 = 1,\\quad \\phi_1 = x,\\quad \\phi_2 = x^2,\\quad \\phi_3 = x^3.\n",
"$$\n",
"Remember that on the interval $[0,3]$ we have the *inner product*\n",
"$$\n",
"(p,q) = \\int_{0}^3p(x)\\,q(x)\\,\\mathrm{d}x\n",
"$$\n",
"and the *norm*\n",
"$$\n",
"\\|p\\| = \\left(\\int_{0}^3 p(x)^2 \\, \\mathrm{d}x\\right)^{1/2}.\n",
"$$\n",
"\n",
"We can now construct orthogonal polynomials by using Gram-Schmidt orthogonalization.\n",
"\n",
"\\begin{align*}\n",
"p_k = \\phi_k - \\sum_{j=0}^{k-1}\\frac{(\\phi_k,p_j)}{\\|p_j\\|^2}\\, p_j\n",
"\\end{align*}\n",
"\n",
"We start out by setting $p_0 = 1$. In order to calculate $p_1$ we first need to calculate\n",
"$$\n",
"(\\phi_1,p_0) = \\int_{0}^3\\phi_1(x)\\,p_0(x)\\, \\mathrm{d}x = \\int_{0}^3 x \\cdot 1 \\, \\mathrm{d}x = \\left[\\frac{x^2}{2}\\right]^3_{0} = \\frac{9}{2}.\n",
"$$\n",
"We also need to calculate\n",
"$$\n",
"\\|p_0\\|^2 = \\int_{0}^3p_0(x)^2 \\, \\mathrm{d}x = \\int_{0}^3 1\\cdot 1 \\mathrm{d}x = 3.\n",
"$$\n",
"Therefore, \n",
"$$\n",
"p_1 = \\phi_1 - \\frac{(\\phi_1,p_0)}{\\|p_0\\|^2}\\,p_0 = \\phi_1 - \\frac{9/2}{3}\\cdot p_0 = \\phi_1 - \\frac{3}{2}p_0 = x-\\frac{3}{2}.\n",
"$$\n",
"\n",
"**a)**\n",
"\n",
"Use Gram-Schmidt orthogonalization to construct $p_2$ and $p_3$."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We can use the Python package SymPy to check our calculations. The code below helps you by defining the inner product and shows how to define polynomials."
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"9/2\n"
]
}
],
"source": [
"from sympy.abc import x\n",
"from sympy import integrate\n",
"\n",
"a=0\n",
"b=3\n",
"\n",
"#Define the inner poduct\n",
"def scp(p,q):\n",
" return integrate(p*q, (x, a, b))\n",
"\n",
"#Define polynomials\n",
"p0 = 1\n",
"p1 = x\n",
"\n",
"#Calculate the inner product and print it.\n",
"print(scp(p0,p1))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"**b)**\n",
"Use the function $\\texttt{scp}$ to check whether the polynomials you calculated\n",
"are in fact **orthogonal**."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"**c)**\n",
"Find the 3 roots of $p_3$.\n",
"\n",
"*Hints:*\n",
"\n",
" - Analytical approach: If you want to do it analytical, use the fact that one root is $$x=\\frac{3}{2}$$\n",
"to find a second order polynomial $\\widetilde{p}_2$ such that $\\widetilde{p}_2(x) \\cdot (x-3/2) = p_3(x)$.\n",
" - Computational approach: \n",
"If you want to use a computational algebra system/symbolic calculator\n",
"you Import the `solve` from `sympy` (Have a look at the\n",
"[solve](https://docs.sympy.org/latest/modules/solvers/solvers.html) submodules.)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"**d)**\n",
"\n",
"Let's denote the three roots of $p_3$ by $x_1,x_2,x_3.$\n",
"\n",
"Construct the three Lagrange polynomials $L_1, \\, L_2,\\, L_3$ satisfying $L_i(x_j) = \\delta_{ij},$ that is\n",
"\n",
"$$\n",
"L_i(x_j) = \n",
"\\begin{cases}\n",
"1,\\quad i=j,\\\\\n",
"0,\\quad i\\neq j.\n",
"\\end{cases}\n",
"$$\n",
"\n",
"Then calculate the weights \n",
"$$\n",
"w_i = \\int_{0}^3L_i(x)\\, \\mathrm{d}x.\n",
"$$\n",
"\n",
"*Hint: You can use the SymPy function $\\texttt{integrate}$ to check your calculations.*"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"**e)** \n",
"\n",
"Now recheck your calculations as follows.\n",
"The Gauss-Legendre rule for the interval $[-1,1]$ with $3$ quadrature points is given by\n",
"\\begin{align*}\n",
"\\{\\hat{x}_i\\}_{i=0}^2 &= \\left\\{-\\sqrt{\\tfrac{3}{5}}, 0, \\sqrt{\\tfrac{3}{5}} \\right\\}\n",
"\\\\\n",
"\\{\\hat{w}_i\\}_{i=0}^2 &= \\left\\{ \\tfrac{5}{9}, \\tfrac{8}{9}, \\tfrac{5}{9} \\right\\}\n",
"\\end{align*}\n",
"Now transfer this quadrature rule to the interval $[0,3]$ and confirm that you get the same quadrature points and weights you computed in 2a)-2d)."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"**f)**\n",
"\n",
"Finally, write down the quadrature rule on the form\n",
"\n",
"$$\n",
"\\mathrm{GQR}[f](0,3) = \\sum_{j=1}^nw_jf(x_j).\n",
"$$\n",
"\n",
"and check that this Gaussian quadrature rule has degree of exactness equal to $5$.\n",
"\n",
"*Hint: Use the `QR` function from the `SimpleQuadrature.ipynb` notebook.*"
]
}
],
"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.9.1"
}
},
"nbformat": 4,
"nbformat_minor": 4
}