File Exchange

image thumbnail

Simpson's 1/3 and 3/8 rules

version 1.2 (3.11 KB) by

SIMPSON: Simpson's rule for quadratic and cubic numerical integration

60 Downloads

Updated

View License

RES = SIMPSON(Y) computes an approximation of the integral of Y via
   Simpson's 1/3 rule (with unit spacing). Simpson's 1/3 rule uses
   quadratic interpolants for numerical integration. To compute the
    integral for spacing different from one, multiply RES by the spacing
    increment.
   
   For vectors, SIMPSON(Y) is the integral of Y. For matrices, SIMPSON(Y)
   is a row vector with the integral over each column. For N-D
    arrays, SIMPSON(Y) works across the first non-singleton dimension.
   
   RES = SIMPSON(X,Y) computes the integral of Y with respect to X using
   Simpson's 1/3 rule. X and Y must be vectors of the same
   length, or X must be a column vector and Y an array whose first
    non-singleton dimension is length(X). SIMPSON operates along this
   dimension. Note that X must be equally spaced for proper execution of
   the 1/3 and 3/8 rules. If X is not equally spaced, the trapezoid rule
   (MATLAB's TRAPZ) is recommended.
   
   RES = SIMPSON(X,Y,DIM) or SIMPSON(Y,DIM) integrates across dimension
    DIM of Y. The length of X must be the same as size(Y,DIM)).
   
    RES = SIMPSON(X,Y,DIM,RULE) can be used to toggle between Simpson's 1/3
    rule and Simpson's 3/8 rule. Simpson's 3/8 rule uses cubic interpolants
    to accomplish the numerical integration. If the default value for DIM
    is desired, assign an empty matrix.
 
    - RULE options
 
        [DEFAULT] '1/3' Simpson's rule for quadratic interpolants
 
                    '3/8' Simpson's rule for cubic interpolants
 
    Examples:
        % Integrate Y = SIN(X)
        x = 0:0.2:pi;
        y = sin(x);
        a = sum(y)*0.2; % Rectangle rule
        b = trapz(x,y); % Trapezoid rule
        c = simpson(x,y,[],'1/3'); % Simpson's 1/3 rule
        d = simpson(x,y,[],'3/8'); % Simpson's 3/8 rule
        e = cos(x(1))-cos(x(end)); % Actual integral
        fprintf('Rectangle Rule: %.15f\n', a)
        fprintf('Trapezoid Rule: %.15f\n', b)
        fprintf('Simpson''s 1/3 Rule: %.15f\n', c)
        fprintf('Simpson''s 3/8 Rule: %.15f\n', d)
        fprintf('Actual Integral: %.15f\n', e)
 
        % http://math.fullerton.edu/mathews/n2003/simpson38rule/Simpson38RuleMod/Links/Simpson38RuleMod_lnk_2.html
        x1 = linspace(0,2,4);
        x2 = linspace(0,2,7);
        x4 = linspace(0,2,13);
        y = @(x) 2+cos(2*sqrt(x));
        format long
        y1 = y(x1); res1 = simpson(x1,y1,[],'3/8'); disp(res1)
        y2 = y(x2); res2 = simpson(x2,y2,[],'3/8'); disp(res2)
        y4 = y(x4); res4 = simpson(x4,y4,[],'3/8'); disp(res4)
 
    Class support for inputs X, Y:
       float: double, single
 
    See also sum, cumsum, trapz, cumtrapz.

Comments and Ratings (3)

Evan Provan

eric ross

Andrew Davis

Thanks for putting this up, it works great.

In the first example, I suppose the step size in x should be 0.01, or the multiplier in the a = ... line should be 0.2. Also, the strcat() function in the disp lines is redundant, since the concatenation is already performed by the [] brackets. May I suggest:
>> fprintf('Rectangle Rule: %.15f\n', a);

Thanks again for this submission.

Updates

1.2

Edited example

1.1

Made inputs and execution congruent with TRAPZ
Help file formatted to MATLAB standard (incl H1 line)
Execution vectorized to accommodate N-dimensional arrays (similar to TRAPZ)
Error checking step included for unequally spaced X

MATLAB Release
MATLAB 7.7 (R2008b)
Acknowledgements

Inspired by: Simpson's Rule Integration

Download apps, toolboxes, and other File Exchange content using Add-On Explorer in MATLAB.

» Watch video