Tutorial 09 - Quadrature
Contents
Tutorial 09 - Quadrature¶
Numerical integration of functions, rectangular rule, trapezoidal rule, Simpson’s rule, Romberg’s method, Gaussian quadrature, Monte-Carlo integration and random number generators.
import numpy as np
import random as rnd
import matplotlib.pyplot as plt
from scipy import integrate, interpolate, special
Newton-Cotes Formulas¶
Midpoint rule¶
Exercise 09.1: Implement the midpoint rule for approximating the definite integral.
def midpoint_rule(f, a, b, N=100):
"""
Calculates definite integral of 1D function using rectangular rule.
Args:
f (function): A function defined on interval [a, b]
a (float): Left-hand side point of the interval
b (float): Right-hand side point of the interval
N (int): Number of subdivisions of the interval [a, b]
Returns:
float: Definite integral
"""
# add your code here
You may evaluate the correctness of your implementation using the scipy.integrate.quad function.
def f(x):
return np.sin(x)
try:
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], \
decimal=5)
except AssertionError as E:
print(E)
else:
print("The implementation is correct.")
Exercise 09.2: Estimate the following definite integral,
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.
# add your code here
Trapezoidal rule¶
Exercise 09.3: Implement the trapezoidal rule for approximating the definite integral.
def trapezoidal_rule(f, a, b, N=100):
"""
Calculates definite integral of 1D function using trapezoidal rule.
Args:
f (function): A function defined on interval [a, b]
a (float): Left-most point of the interval
b (float): Right-most point of the interval
N (int): Number of subdivisions of the interval [a, b]
Returns:
float: Definite integral
"""
# add your code here
You may evaluate the correctness of your implementation using the scipy.integrate.quad function.
def f(x):
return np.sin(x)
try:
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], \
decimal=5)
except AssertionError as E:
print(E)
else:
print("The implementation is correct.")
Exercise 09.4: Estimate the following definite integral,
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.
# add your code here
Simpson’s quadratic rule¶
Exercise 09.5: Implement the Simpsons’s quadratic rule for approximating the definite integral.
def simpsons_quadratic_rule(f, a, b, N=100):
"""
Calculates definite integral of 1D function using Simpson's 1/3 rule.
Args:
f (function): A function defined on interval [a, b]
a (float): Left-hand side point of the interval
b (float): Right-hand side point of the interval
N (int): Number of subdivisions of the interval [a, b]
Returns:
float: Definite integral
"""
# add your code here
You may evaluate the correctness of your implementation using the scipy.integrate.quad function.
def f(x):
return np.sin(x)
try:
np.testing.assert_almost_equal(simpsons_quadratic_rule(f, 0.0, 2.0 * np.pi, 10000), \
integrate.quad(f, 0.0, 2.0 * np.pi)[0], decimal=5)
except AssertionError as E:
print(E)
else:
print("The implementation is correct.")
Exercise 09.6: Estimate the following definite integral,
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 function). Find a value \( N \) which guarantees that the approximation is within \( 10^{-5} \) of the exact value.
# add your code here
Note: The integral in Exercise 09.6 is called the Fresnel integral.
Simpson’s cubic rule¶
Exercise 09.7: Implement the Simpson’s cubic rule for approximating the definite integral.
def simpsons_cubic_rule(f, a, b, n=100):
"""
Calculates definite integral of 1D function using Simpson's 3/8 rule.
Args:
f (function): A function defined on interval [a, b]
a (float): Left-hand side point of the interval
b (float): Right-hand side point of the interval
n (int): Number of subdivisions of the interval [a, b]
Returns:
float: Definite integral
"""
# add your code here
You may evaluate the correctness of your implementation using the scipy.integrate.quad function.
def f(x):
return np.sin(x)
try:
np.testing.assert_almost_equal(simpsons_cubic_rule(f, 0.0, 2.0 * np.pi, 10000), \
integrate.quad(f, 0.0, 2.0 * np.pi)[0], decimal=5)
except AssertionError as E:
print(E)
else:
print("The implementation is correct.")
Exercise 09.8: Estimate the following definite integral,
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 function). Find a value \( N \) which guarantees that the approximation is within \( 10^{-5} \) of the exact value.
# add your code here
Note: The integral in Exercise 09.8 is called the error function.
Romberg’s method¶
Exercise 09.9: Implement the Romberg’s method for approximating the definite integral.
def rombergs_method(f, a, b, m=10):
"""
Calculates definite integral of 1D function using Romberg's method.
Args:
f (function): A function defined on interval [a, b]
a (float): Left-hand side point of the interval
b (float): Right-hand side point of the interval
m (int): Number of extrapolations
Returns:
float: Definite integral
"""
# add your code here
You may evaluate the correctness of your implementation using the scipy.integrate.quad function.
def f(x):
return np.sin(x)
try:
np.testing.assert_almost_equal(rombergs_method(f, 0.0, 2.0 * np.pi, 10), \
integrate.quad(f, 0.0, 2.0 * np.pi)[0], decimal=5)
except AssertionError as E:
print(E)
else:
print("The implementation is correct.")
Gaussian quadrature¶
Exercise 09.10: Write a function that calculates a Legendre polynomial of degree \( n \) using the Bonnet’s recursion formula,
def legendre_polynomial(n):
"""
Calculates a Legendre polynomial of degree n using recursive formula.
Args:
n (int): Degree of the polynomial
Returns:
numpy.poly1d: The Legendre polynomial of degree n
"""
# add your code here
You may evaluate the correctness of your implementation using the scipy.special.legendre function.
n = 10
x = np.linspace(-1.0, 1.0, 1000)
try:
np.testing.assert_array_almost_equal(legendre_polynomial(n)(x), special.legendre(n)(x), decimal=7)
except AssertionError as E:
print(E)
else:
print("The implementation is correct.")
Exercise 09.11: Plot the first 5 Legendre polynomials on \( x \in [-1, 1] \).
# add your code here
Gauss-Legendre quadrature rule¶
Gauss-Legendre quadrature rule
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 function.
def gauss_legendre_quadrature(f, a, b, N=3):
"""
Calculates definite integral of 1D function using Gaussian quadrature rule.
Args:
f (function): A function defined on interval [a, b]
a (float): Left-hand side point of the interval
b (float): Right-hand side point of the interval
N (int): Degree of used polynomial
Returns:
float: Definite integral
"""
# add your code here
You may evaluate the correctness of your implementation using the scipy.integrate.quad function.
def f(x):
return np.sin(x)
try:
np.testing.assert_almost_equal(gauss_legendre_quadrature(f, 0.0, 2.0 * np.pi, 3), \
integrate.quad(f, 0.0, 2.0 * np.pi)[0], decimal=5)
except AssertionError as E:
print(E)
else:
print("The implementation is correct.")
Monte Carlo integration¶
Exercise 09.13: Implement the Monte Carlo method for approximating the definite integral.
def monte_carlo_integration(f, a, b, c, d, N=500):
"""
Calculates definite integral of 1D function using Monte Carlo method.
Args:
f (function): A function defined on interval [a, b]
a (float): Left-hand side point of the interval
b (float): Right-hand side point of the interval
N (int): Number of random points generated
Returns:
float: Definite integral
"""
# add your code here
Exercise 09.14: Estimate the following definite integral,
using the Monte Carlo method with \( N = 1000 \) samples. Calculate the absolute error of the approximation (the exact value of the integral is \( 0 \)).
# add your code here