{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Tutorial 09 - Quadrature"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Numerical integration of functions, rectangular rule, trapezoidal rule, Simpson's rule, Romberg's method, Gaussian quadrature, Monte-Carlo integration and random number generators. "
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"import numpy as np\n",
"import random as rnd\n",
"import matplotlib.pyplot as plt\n",
"from scipy import integrate, interpolate, special"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Newton-Cotes Formulas"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Midpoint rule"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"[Midpoint rule](https://en.wikipedia.org/wiki/Riemann_sum#Midpoint_rule)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"\n",
"**Exercise 09.1:** Implement the midpoint rule for approximating the definite integral.\n",
" \n",
"
"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"def midpoint_rule(f, a, b, N=100):\n",
" \"\"\"\n",
" Calculates definite integral of 1D function using rectangular rule.\n",
" Args:\n",
" f (function): A function defined on interval [a, b]\n",
" a (float): Left-hand side point of the interval\n",
" b (float): Right-hand side point of the interval\n",
" N (int): Number of subdivisions of the interval [a, b]\n",
" Returns:\n",
" float: Definite integral\n",
" \"\"\"\n",
"\n",
" # add your code here\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"You may evaluate the correctness of your implementation using the [`scipy.integrate.quad`](https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.quad.html) function."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"def f(x):\n",
" return np.sin(x)\n",
"try:\n",
" np.testing.assert_almost_equal(midpoint_rule(f, 0.0, 2.0 * np.pi, 10000), integrate.quad(f, 0.0, 2.0 * np.pi)[0], \\\n",
" decimal=5)\n",
"except AssertionError as E:\n",
" print(E)\n",
"else:\n",
" print(\"The implementation is correct.\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"\n",
"**Exercise 09.2:** Estimate the following definite integral,\n",
" \n",
"$$\n",
"\\int_0^1 \\frac{4}{1 + x^2} dx\n",
"$$\n",
"\n",
"using the midpoint rule with a partition of size $ N = 10 $. Calculate the absolute error of the approximation (the exact value of the integral is $ \\pi $). Find a value $ N $ which guarantees that the approximation is within $ 10^{-5} $ of the exact value.\n",
" \n",
"
"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"\n",
"# add your code here\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Trapezoidal rule"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"[Trapezoidal rule](https://en.wikipedia.org/wiki/Trapezoidal_rule)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"\n",
"**Exercise 09.3:** Implement the trapezoidal rule for approximating the definite integral.\n",
" \n",
"
"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"def trapezoidal_rule(f, a, b, N=100):\n",
" \"\"\"\n",
" Calculates definite integral of 1D function using trapezoidal rule.\n",
" Args:\n",
" f (function): A function defined on interval [a, b]\n",
" a (float): Left-most point of the interval\n",
" b (float): Right-most point of the interval\n",
" N (int): Number of subdivisions of the interval [a, b]\n",
" Returns:\n",
" float: Definite integral\n",
" \"\"\"\n",
" \n",
" # add your code here\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"You may evaluate the correctness of your implementation using the [`scipy.integrate.quad`](https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.quad.html) function."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"def f(x):\n",
" return np.sin(x)\n",
"try:\n",
" np.testing.assert_almost_equal(trapezoidal_rule(f, 0.0, 2.0 * np.pi, 10000), integrate.quad(f, 0.0, 2.0 * np.pi)[0], \\\n",
" decimal=5)\n",
"except AssertionError as E:\n",
" print(E)\n",
"else:\n",
" print(\"The implementation is correct.\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"\n",
"**Exercise 09.4:** Estimate the following definite integral,\n",
" \n",
"$$\n",
"\\int_1^2 \\frac{1}{x} dx\n",
"$$\n",
"\n",
"using the trapezoidal rule with a partition of size $ N = 10 $. Calculate the absolute error of the approximation (the exact value of the integral is $ \\ln 2 $). Find a value $ N $ which guarantees that the approximation is within $ 10^{-5} $ of the exact value.\n",
" \n",
"
"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"\n",
"# add your code here\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Simpson's quadratic rule"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"[Simpson's quadratic rule](https://en.wikipedia.org/wiki/Simpson%27s_rule#Simpson's_1/3_rule)\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"\n",
"**Exercise 09.5:** Implement the Simpsons's quadratic rule for approximating the definite integral.\n",
" \n",
"
"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"def simpsons_quadratic_rule(f, a, b, N=100):\n",
" \"\"\"\n",
" Calculates definite integral of 1D function using Simpson's 1/3 rule.\n",
" Args:\n",
" f (function): A function defined on interval [a, b]\n",
" a (float): Left-hand side point of the interval\n",
" b (float): Right-hand side point of the interval\n",
" N (int): Number of subdivisions of the interval [a, b]\n",
" Returns:\n",
" float: Definite integral\n",
" \"\"\"\n",
" \n",
" # add your code here\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"You may evaluate the correctness of your implementation using the [`scipy.integrate.quad`](https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.quad.html) function."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"def f(x):\n",
" return np.sin(x)\n",
"try:\n",
" np.testing.assert_almost_equal(simpsons_quadratic_rule(f, 0.0, 2.0 * np.pi, 10000), \\\n",
" integrate.quad(f, 0.0, 2.0 * np.pi)[0], decimal=5)\n",
"except AssertionError as E:\n",
" print(E)\n",
"else:\n",
" print(\"The implementation is correct.\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"\n",
"**Exercise 09.6:** Estimate the following definite integral,\n",
" \n",
"$$\n",
"\\int_0^2 \\sin \\left( \\frac{\\pi x^2}{2} \\right) dx\n",
"$$\n",
"\n",
"using the Simpson's quadratic rule with a partition of size $ N = 10 $. Calculate the absolute error of the approximation (the exact value of the integral can be obtained by the [`scipy.special.fresnel`](https://docs.scipy.org/doc/scipy/reference/generated/scipy.special.fresnel.html) function). Find a value $ N $ which guarantees that the approximation is within $ 10^{-5} $ of the exact value.\n",
" \n",
"
"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
" \n",
"# add your code here\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
" \n",
"**Note:** The integral in Exercise 09.6 is called the [Fresnel integral](https://en.wikipedia.org/wiki/Fresnel_integral).\n",
" \n",
"
"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Simpson's cubic rule"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"[Simpson's cubic rule](https://en.wikipedia.org/wiki/Simpson%27s_rule#Simpson's_3/8_rule)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"\n",
"**Exercise 09.7:** Implement the Simpson's cubic rule for approximating the definite integral.\n",
" \n",
"
"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"def simpsons_cubic_rule(f, a, b, n=100):\n",
" \"\"\"\n",
" Calculates definite integral of 1D function using Simpson's 3/8 rule.\n",
" Args:\n",
" f (function): A function defined on interval [a, b]\n",
" a (float): Left-hand side point of the interval\n",
" b (float): Right-hand side point of the interval\n",
" n (int): Number of subdivisions of the interval [a, b]\n",
" Returns:\n",
" float: Definite integral\n",
" \"\"\"\n",
" \n",
" # add your code here\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"You may evaluate the correctness of your implementation using the [`scipy.integrate.quad`](https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.quad.html) function."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"def f(x):\n",
" return np.sin(x)\n",
"try:\n",
" np.testing.assert_almost_equal(simpsons_cubic_rule(f, 0.0, 2.0 * np.pi, 10000), \\\n",
" integrate.quad(f, 0.0, 2.0 * np.pi)[0], decimal=5)\n",
"except AssertionError as E:\n",
" print(E)\n",
"else:\n",
" print(\"The implementation is correct.\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"\n",
"**Exercise 09.8:** Estimate the following definite integral,\n",
" \n",
"$$\n",
"\\int_0^1 \\frac{2}{\\sqrt{\\pi}} e^{-x^2} dx\n",
"$$\n",
"\n",
"using the Simpson's cubic rule with a partition of size $ N = 10 $. Calculate the absolute error of the approximation (the exact value of the integral can be obtained by the [`scipy.special.erf`](https://docs.scipy.org/doc/scipy/reference/generated/scipy.special.erf.html) function). Find a value $ N $ which guarantees that the approximation is within $ 10^{-5} $ of the exact value.\n",
" \n",
" \n",
"
"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"\n",
"# add your code here\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
" \n",
"**Note:** The integral in Exercise 09.8 is called the [error function](https://en.wikipedia.org/wiki/Error_function).\n",
" \n",
"
"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Romberg's method"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"[Romberg's method](https://en.wikipedia.org/wiki/Romberg%27s_method)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"\n",
"**Exercise 09.9:** Implement the Romberg's method for approximating the definite integral.\n",
" \n",
"
"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"def rombergs_method(f, a, b, m=10):\n",
" \"\"\"\n",
" Calculates definite integral of 1D function using Romberg's method.\n",
" Args:\n",
" f (function): A function defined on interval [a, b]\n",
" a (float): Left-hand side point of the interval\n",
" b (float): Right-hand side point of the interval\n",
" m (int): Number of extrapolations\n",
" Returns:\n",
" float: Definite integral\n",
" \"\"\"\n",
" \n",
" # add your code here\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"You may evaluate the correctness of your implementation using the [`scipy.integrate.quad`](https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.quad.html) function."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"def f(x):\n",
" return np.sin(x)\n",
"try:\n",
" np.testing.assert_almost_equal(rombergs_method(f, 0.0, 2.0 * np.pi, 10), \\\n",
" integrate.quad(f, 0.0, 2.0 * np.pi)[0], decimal=5)\n",
"except AssertionError as E:\n",
" print(E)\n",
"else:\n",
" print(\"The implementation is correct.\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Gaussian quadrature"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"\n",
"**Exercise 09.10:** Write a function that calculates a [Legendre polynomial](https://en.wikipedia.org/wiki/Legendre_polynomials) of degree $ n $ using the Bonnet’s recursion formula,\n",
" \n",
"$$\n",
"P_n(x) = \\frac{2n - 1}{n} x P_{n-1}(x) - \\frac{n-1}{n} P_{n-2}(x), \\quad P_0(x) = 1, \\ P_1(x) = x.\n",
"$$\n",
" \n",
"
"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"def legendre_polynomial(n):\n",
" \"\"\"\n",
" Calculates a Legendre polynomial of degree n using recursive formula.\n",
" Args:\n",
" n (int): Degree of the polynomial\n",
" Returns:\n",
" numpy.poly1d: The Legendre polynomial of degree n\n",
" \"\"\"\n",
"\n",
" # add your code here\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"You may evaluate the correctness of your implementation using the [`scipy.special.legendre`](https://docs.scipy.org/doc/scipy/reference/generated/scipy.special.legendre.html) function."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"n = 10\n",
"x = np.linspace(-1.0, 1.0, 1000)\n",
"try:\n",
" np.testing.assert_array_almost_equal(legendre_polynomial(n)(x), special.legendre(n)(x), decimal=7)\n",
"except AssertionError as E:\n",
" print(E)\n",
"else:\n",
" print(\"The implementation is correct.\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"\n",
"**Exercise 09.11:** Plot the first 5 Legendre polynomials on $ x \\in [-1, 1] $.\n",
" \n",
"
\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
" \n",
"# add your code here\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Gauss-Legendre quadrature rule"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"[Gauss-Legendre quadrature rule](https://en.wikipedia.org/wiki/Gaussian_quadrature#Gauss%E2%80%93Legendre_quadrature)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"\n",
"**Exercise 09.12:** Implement the Gauss–Legendre quadrature rule for approximating the definite integral. You may calculate the roots of Legendre polynomials and their weights using the [`numpy.polynomial.legendre.leggauss`](https://numpy.org/doc/stable/reference/generated/numpy.polynomial.legendre.leggauss.html) function. \n",
" \n",
"
"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"def gauss_legendre_quadrature(f, a, b, N=3):\n",
" \"\"\"\n",
" Calculates definite integral of 1D function using Gaussian quadrature rule.\n",
" Args:\n",
" f (function): A function defined on interval [a, b]\n",
" a (float): Left-hand side point of the interval\n",
" b (float): Right-hand side point of the interval\n",
" N (int): Degree of used polynomial\n",
" Returns:\n",
" float: Definite integral\n",
" \"\"\"\n",
" \n",
" # add your code here\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"You may evaluate the correctness of your implementation using the [`scipy.integrate.quad`](https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.quad.html) function."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"def f(x):\n",
" return np.sin(x)\n",
"try:\n",
" np.testing.assert_almost_equal(gauss_legendre_quadrature(f, 0.0, 2.0 * np.pi, 3), \\\n",
" integrate.quad(f, 0.0, 2.0 * np.pi)[0], decimal=5)\n",
"except AssertionError as E:\n",
" print(E)\n",
"else:\n",
" print(\"The implementation is correct.\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Monte Carlo integration"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"[Monte Carlo integration](https://en.wikipedia.org/wiki/Monte_Carlo_integration)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"\n",
"**Exercise 09.13:** Implement the Monte Carlo method for approximating the definite integral. \n",
" \n",
"
"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"def monte_carlo_integration(f, a, b, c, d, N=500):\n",
" \"\"\"\n",
" Calculates definite integral of 1D function using Monte Carlo method.\n",
" Args:\n",
" f (function): A function defined on interval [a, b]\n",
" a (float): Left-hand side point of the interval\n",
" b (float): Right-hand side point of the interval\n",
" N (int): Number of random points generated\n",
" Returns:\n",
" float: Definite integral\n",
" \"\"\"\n",
" \n",
" # add your code here\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"\n",
"**Exercise 09.14:** Estimate the following definite integral,\n",
" \n",
"$$\n",
"\\int_0^{4 \\pi} \\sin(x) dx\n",
"$$\n",
"\n",
"using the Monte Carlo method with $ N = 1000 $ samples. Calculate the absolute error of the approximation (the exact value of the integral is $ 0 $).\n",
" \n",
" \n",
"
"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"scrolled": true
},
"outputs": [],
"source": [
"\n",
"# add your code here\n"
]
}
],
"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.8.5"
}
},
"nbformat": 4,
"nbformat_minor": 2
}