diff options
Diffstat (limited to 'tutorials/module_3/4_numerical_integration.md')
| -rw-r--r-- | tutorials/module_3/4_numerical_integration.md | 177 |
1 files changed, 162 insertions, 15 deletions
diff --git a/tutorials/module_3/4_numerical_integration.md b/tutorials/module_3/4_numerical_integration.md index 6e534b5..78b0328 100644 --- a/tutorials/module_3/4_numerical_integration.md +++ b/tutorials/module_3/4_numerical_integration.md @@ -47,40 +47,187 @@ $$ It is simple to implement and works especially well for tabulated data. Like the midpoint rule, its accuracy is of order $O(h^2)$. + + +```python +import numpy as np + +a = 0 +b = np.pi +n = 11 +h = (b - a) / (n - 1) +x = np.linspace(a, b, n) +f = np.sin(x) + +I_trap = (h/2)*(f[0] + 2 * sum(f[1:n-1]) + f[n-1]) +err_trap = 2 - I_trap + +print(I_trap) +print(err_trap) +``` + + ### Simpson’s Rule Simpson’s rules use polynomial interpolation to achieve higher accuracy. -- **Simpson’s 1/3 Rule (order $O(h^4)$)** +**Simpson’s 1/3 Rule (order $O(h^4)$)** Requires an even number of sub-intervals $n$: $$ I \approx \frac{h}{3}\Big[f(x_0) + 4\sum_{\text{odd } i} f(x_i) + 2\sum_{\text{even } i<n} f(x_i) + f(x_n)\Big]. $$ +```python +import numpy as np + +a = 0 +b = np.pi +n = 11 +h = (b - a) / (n - 1) +x = np.linspace(a, b, n) +f = np.sin(x) + +I_simp = (h/3) * (f[0] + 2*sum(f[:n-2:2]) \ + + 4*sum(f[1:n-1:2]) + f[n-1]) +err_simp = 2 - I_simp -- **Simpson’s 3/8 Rule** - Based on cubic interpolation across three subintervals: +print(I_simp) +print(err_simp) +``` + + +**Simpson’s 3/8 Rule** +Based on cubic interpolation across three subintervals: $$ I \approx \frac{3h}{8}\Big[f(x_0) + 3f(x_1) + 3f(x_2) + f(x_3)\Big]. $$ - Often used to handle the final three panels when \(n\) is not divisible by two. +Often used to handle the final three panels when \(n\) is not divisible by two. + +```python +import numpy as np + +a = 0 +b = np.pi +n = 11 +h = (b - a) / (n - 1) +x = np.linspace(a, b, n) +f = np.sin(x) + +I_simp = (h/3) * (f[0] + 2*sum(f[:n-2:2]) \ + + 4*sum(f[1:n-1:2]) + f[n-1]) +err_simp = 2 - I_simp + +print(I_simp) +print(err_simp) +``` ### Romberg Integration Romberg integration refines trapezoidal approximations using Richardson extrapolation to cancel error terms. This results in very rapid convergence for smooth functions, making it a powerful method once the trapezoidal rule is implemented. +```python +import numpy as np + +a = 0.0 +b = 1.0 + +f = lambda x: np.exp(-x**2) + +iter_ = 0 +iter_max = 50 +err = 1.0 +tol = 1e-14 + +I = np.zeros((iter_max + 2, iter_max + 2), dtype=float) + +def trap_func(func, a, b, n): + """Composite trapezoidal rule with n subintervals on [a, b].""" + h = (b - a) / n + x = np.linspace(a, b, n + 1) + y = func(x) + return h * (0.5 * y[0] + y[1:-1].sum() + 0.5 * y[-1]) + + +n = 1 +I[0, 0] = trap_func(f, a, b, n) + +# ---- Romberg loop ---- +while (iter_ < iter_max) and (err > tol): + iter_ += 1 + n = 2 ** iter_ + I[iter_, 0] = trap_func(f, a, b, n) + for k in range(1, iter_ + 1): + j = iter_ - k + I[j, k] = (4**k * I[j + 1, k - 1] - I[j, k - 1]) / (4**k - 1) + if iter_ >= 1: + err = abs((I[0, iter_] - I_]() +``` + ### Gaussian Quadrature -Gaussian quadrature chooses evaluation points and weights optimally (based on Legendre -polynomials) to maximize accuracy. With $n$ evaluation points, Gaussian quadrature is exact for polynomials of degree up to $2n-1$. This makes it extremely efficient for smooth integrands. +Gaussian quadrature chooses evaluation points and weights optimally (based on Legendre polynomials) to maximize accuracy. With $n$ evaluation points, Gaussian quadrature is exact for polynomials of degree up to $2n-1$. This makes it extremely efficient for smooth integrands. + +```python +import numpy as np + +a, b = 0.0, 1.0 # integration limits [a, b] +n_list = [2, 3, 4, 5, 8, 16] # quadrature orders to try + +for n in n_list: + xi, wi = np.polynomial.legendre.leggauss(n) + x = 0.5*(b - a)*xi + 0.5*(b + a) + f_vals = np.exp(-x**2) # f(x) = e^(-x^2) + + I_approx = 0.5*(b - a) * np.sum(wi * f_vals) + print(f"{n:<3}| {I_approx:.12f}") + +``` + + + +```python +#!/usr/bin/env python3 +# Gauss Quadrature using Gauss-Legendre method (manual 2-point & 3-point) + +import numpy as np +from scipy.integrate import quad # for reference "exact" integral + +# ---- Establish a, b, and the function ---- +a = 0.0 +b = 1.0 +f = lambda x: np.exp(-x**2) # example: f(x) = e^(-x^2) + +# ---- Common factor ---- +dx = (b - a) / 2.0 + +# ---- Two-point Gauss-Legendre ---- +x0 = ((b + a) + (b - a) * (-1/np.sqrt(3))) / 2 +x1 = ((b + a) + (b - a) * ( 1/np.sqrt(3))) / 2 +I_2Point = (f(x0) + f(x1)) * dx +print(f"2-point Gauss-Legendre: {I_2Point:.10f}") + +# ---- Three-point Gauss-Legendre ---- +c0, c1, c2 = 5/9, 8/9, 5/9 +x0 = ((b + a) + (b - a) * (-np.sqrt(3/5))) / 2 +x1 = ((b + a) + (b - a) * (0.0)) / 2 +x2 = ((b + a) + (b - a) * ( np.sqrt(3/5))) / 2 +I_3Point = (c0*f(x0) + c1*f(x1) + c2*f(x2)) * dx +print(f"3-point Gauss-Legendre: {I_3Point:.10f}") + +# ---- Actual value using scipy.integrate.quad ---- +I_actual, _ = quad(f, a, b) +print(f"Actual (scipy.quad): {I_actual:.10f}") + +# ---- Errors ---- +err2 = abs(I_actual - I_2Point) +err3 = abs(I_actual - I_3Point) +print(f"Error (2-point): {err2:.2e}") +print(f"Error (3-point): {err3:.2e}") + +``` + + ---- -## Taking this further -- **Romberg Integration:** Demonstrates how extrapolation accelerates convergence of trapezoidal approximations. -- **Gaussian Quadrature:** Introduces optimal integration points and highlights efficiency for -polynomial and smooth functions. -- **Simpson’s Rules:** Show how higher-order polynomial interpolation improves accuracy. -These methods will be implemented and compared in subsequent assignments to build a deeper understanding of numerical integration accuracy and efficiency. --- @@ -121,12 +268,12 @@ Students should compare results for increasing $n$ and observe how the error dec ## Gaussian Quadrature -Implement a Python function for two-point and three-point Gauss–Legendre quadrature over an arbitrary interval $[a,b]$. Verify exactness for polynomials up to the appropriate degree and compare performance against the trapezoidal rule on oscillatory test functions. +Write a Python function for two-point and three-point Gauss–Legendre quadrature over an arbitrary interval $[a,b]$. Verify exactness for polynomials up to the appropriate degree and compare performance against the trapezoidal rule on oscillatory test functions. --- ## Simpson’s 1/3 Rule -Implement the composite Simpson’s 1/3 rule. Test its accuracy on smooth functions and compare its performance to the trapezoidal rule and Gaussian quadrature. Document error trends and discuss cases where Simpson’s method is preferable. +Code the composite Simpson’s 1/3 rule. Test its accuracy on smooth functions and compare its performance to the trapezoidal rule and Gaussian quadrature. Document error trends and discuss cases where Simpson’s method is preferable. --- |