将梯形法改写为辛普森法的MATLAB实现及积分程序开发问询

Got it, let's adapt your trapezoidal code to implement the composite Simpson's rule! The key idea here is leveraging the relationship between trapezoidal and Simpson's integrals to avoid redundant calculations, which keeps the code efficient and aligned with your existing structure.

Key Background #

Composite Simpson's rule approximates the integral using quadratic polynomials over pairs of subintervals (so the total number of subintervals must be even). A handy recursive formula lets us compute Simpson's integral from trapezoidal estimates:
S_n = (4*T_n - T_{n/2})/3
where T_n is the trapezoidal estimate with n subintervals, and T_{n/2} is the estimate from the previous (coarser) step.

Rewritten Code #

Here's the modified function for composite Simpson's rule, keeping your original function signature and output structure:

function [int, flag, stats] = simpComp(f, a, b, tol, hMin)
    % Initialise variables
    h = b - a;
    n = 1;
    % Compute initial trapezoidal estimates for the first Simpson's step
    T_prev = h / 2 * (f(a) + f(b)); % T_1 (1 subinterval)
    h = h / 2;
    n = 2 * n;
    x_mid = a + h;
    T_curr = T_prev / 2 + h * f(x_mid); % T_2 (2 subintervals)
    int = (4*T_curr - T_prev)/3; % Initial Simpson's estimate S_1
    flag = 1;
    
    % Initialise stats if requested
    if nargout == 3
        stats = struct('totalErEst', [], 'totalNrIntervals', [], 'nodesList', []);
        stats.totalNrIntervals = [stats.totalNrIntervals, n];
        stats.nodesList = [stats.nodesList, a:h:b];
        % Estimate error (Simpson's error is ~ (S_curr - S_prev)/15 for next step)
        stats.totalErEst = [stats.totalErEst, abs(int - (4*T_curr - T_prev)/3)]; % Initial error placeholder
    end
    
    while h > hMin
        T_prev = T_curr;
        h = h / 2;
        n = 2 * n;
        
        if h < eps % Check if h is effectively zero to avoid infinite loop
            break;
        end
        
        % Compute new trapezoidal estimate by adding midpoints of new subintervals
        x_new = a + h : 2*h : b - h;
        T_curr = T_prev / 2 + h * sum(f(x_new));
        
        % Update Simpson's integral using the recursive formula
        int_prev = int;
        int = (4*T_curr - T_prev)/3;
        
        % Update flag if tolerance is met
        if abs(int - int_prev) < tol
            flag = 0;
            if nargout == 3
                stats.totalErEst = [stats.totalErEst, abs(int - int_prev)];
                stats.totalNrIntervals = [stats.totalNrIntervals, n];
                stats.nodesList = [stats.nodesList, a:h:b];
            end
            break;
        end
        
        % Update stats if requested
        if nargout == 3
            stats.totalErEst = [stats.totalErEst, abs(int - int_prev)];
            stats.totalNrIntervals = [stats.totalNrIntervals, n];
            stats.nodesList = [stats.nodesList, a:h:b];
        end
    end
    
    % If loop exited due to hMin, set flag to indicate no convergence to tol
    if flag == 1 && h <= hMin
        flag = 2;
    end
end

Key Changes from Trapezoidal Code #

• Initialization: We start by computing the first two trapezoidal estimates (T_1 and T_2) to get the initial Simpson's integral using the recursive formula.
• Recursive Update: Instead of updating the trapezoidal integral directly, we compute new trapezoidal estimates and convert them to Simpson's integrals using (4*T_curr - T_prev)/3.
• Error Checking: We track the difference between consecutive Simpson's estimates to check against the tolerance tol.
• Stats Tracking: Updated to store Simpson-specific interval counts, node lists, and error estimates.

Notes #

• Simpson's rule converges faster (O(h⁴)) than the trapezoidal rule (O(h²)), so you'll likely meet the tolerance tol with fewer intervals.
• The code ensures the number of subintervals n is always even, which is required for composite Simpson's rule.

内容的提问来源于stack exchange，提问作者M. Hartog