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
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,
$$ \int_0^1 \frac{4}{1 + x^2} dx $$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
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,
$$ \int_1^2 \frac{1}{x} dx $$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
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,
$$ \int_0^2 \sin \left( \frac{\pi x^2}{2} \right) dx $$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.
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,
$$ \int_0^1 \frac{2}{\sqrt{\pi}} e^{-x^2} dx $$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.
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.")
Exercise 09.10: Write a function that calculates a Legendre polynomial of degree $ n $ using the Bonnet’s recursion formula,
$$ 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. $$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
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.")
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,
$$ \int_0^{4 \pi} \sin(x) dx $$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