Advanced Numerical Integration Methods: Gaussian Quadrature Implementation in MATLAB

# 2.1 Derivation of Gaussian Quadrature Formula The Gaussian quadrature formula is derived from the theory of orthogonal polynomials. For a function `f(x)` on the interval `[a, b]`, its sequence of orthogonal polynomials `{p_n(x)}` satisfies: ``` ∫[a, b] p_i(x) p_j(x) w(x) dx = δ_{i, j} ``` where `w(x)` is the weight function on the interval `[a, b]`, and `δ_{i, j}` is the Kronecker delta. For a polynomial `p_n(x)` of degree `n` on the interval `[a, b]`, the Gaussian quadrature formula is: ``` ∫[a, b] f(x) dx ≈ ∑_{i=1}^n w_i f(x_i) ``` where `x_i` are the `n` orthogonal roots of `p_n(x)`, and `w_i` are the corresponding weights. The derivation process of the Gaussian quadrature formula is as follows: 1. Construct the sequence of orthogonal polynomials `{p_n(x)}`. 2. Find the orthogonal roots `{x_i}` of the orthogonal polynomial `p_n(x)`. 3. Determine the corresponding weights `{w_i}` for the orthogonal roots. 4. Substitute the values of `f(x)` at the orthogonal roots into the Gaussian quadrature formula to obtain an approximate integral value. # 2.1 Derivation of Gaussian Quadrature Formula The Gaussian quadrature formula is derived from the theory of orthogonal polynomials. For a function `f(x)` on the interval `[a, b]`, its Gaussian quadrature formula is: ``` ∫[a, b] f(x) dx ≈ ∑[i=1, n] w[i] f(x[i]) ``` where: * `n` is the number of quadrature points * `w[i]` are the weight coefficients * `x[i]` are the integration nodes **Derivation Process:** 1. **Construct Orthogonal Polynomials:** For the interval `[a, b]`, construct a set of orthogonal polynomials `p[i](x)` that satisfy: ``` ∫[a, b] p[i](x) p[j](x) dx = δ[i, j] ``` where `δ[i, j]` is the Kronecker delta. 2. **Construct the Quadrature Formula:** Expand the function `f(x)` over the interval `[a, b]` using orthogonal polynomials: ``` f(x) = ∑[i=0, ∞] c[i] p[i](x) ``` where `c[i]` are the expansion coefficients. Substitute the expansion into the integral formula to get: ``` ∫[a, b] f(x) dx = ∫[a, b] ∑[i=0, ∞] c[i] p[i](x) dx ``` Due to the property of orthogonal polynomials, only a finite number of terms are non-zero, thus the expansion can be truncated: ``` ∫[a, b] f(x) dx ≈ ∫[a, b] ∑[i=0, n-1] c[i] p[i](x) dx ``` where `n` is the truncation order. 3. **Determine Integration Nodes and Weight Coefficients:** Select `n` integration nodes `x[i]` such that the following system of equations holds: ``` ∫[a, b] p[j](x) dx = ∑[i=1, n] w[i] p[j](x[i]) ``` Solving the system of equations gives the integration nodes and weight coefficients. **Derivation Results:** For the interval `[a, b]`, the integration nodes and weight coefficients of the Gaussian quadrature formula are: ``` x[i] = (b - a) / 2 * cos(π (i - 1/2) / n) + (a + b) / 2 w[i] = (b - a) / 2 * sin(π (i - 1/2) / n) / sin(π (i - 1/2) / 2n) ``` # 3.1 Gaussian Quadrature Functions in MATLAB MATLAB provides an inbuilt function `quadgk` for computing one-dimensional integrals. This function uses an adaptive Gaussian quadrature algorithm, which automatically selects the quadrature points and weights to achieve high precision. ``` % Syntax integral = quadgk(@fun, a, b, tol) % Parameter Description - fun: Handle to the integrand function - a: Lower limit of integration - b: Upper limit of integration - tol: Tolerance, specifies the accuracy of integration ``` **Code Logic Analysis:** 1. The `quadgk` function accepts four parameters: the handle to the integrand function, lower and upper limits of integration, and tolerance. 2. The function internally uses an adaptive Gaussian quadrature algorithm, which automatically selects the quadrature points and weights based on the complexity of the integrand and the tolerance requirements. 3. The algorithm employs a divide-and-conquer strategy, subdividing the integration interval until the tolerance requirements are met. 4. The function returns the result of the integral, representing the value of the integrand over the given interval. **Example:** ``` % Integrand function fun = @(x) exp(-x.^2); % Integration interval a = -1; b = 1; % Tolerance tol = 1e-6; % Compute integral integral = quadgk(fun, a, b, tol); % Output result fprintf('Integral result: %f\n', integral); ``` **Output:** ``` Integral result: 1.*** ``` ## 3.2 Parameter Selection for Gaussian Quadrature Method The accuracy of the Gaussian quadrature method is affected by the choice of quadrature points and weights. The `quadgk` function in MATLAB uses an adaptive algorithm, automatically selecting parameters based on the complexity of the integrand and the tolerance requirements. **Adaptive Algorithm:** The adaptive algorithm is based on the following principles: 1. Divide the integration interval into small*** *** ***pare the integral values of adjacent subintervals, and if the difference exceeds the tolerance, further subdivide the subinterval. 4. Repe