File Exchange

image thumbnail

Simpson's 1/3 and 3/8 rules

version 1.2.0.0 (3.11 KB) by Jered Wells
SIMPSON: Simpson's rule for quadratic and cubic numerical integration

25 Downloads

Updated 23 Mar 2012

View Version History

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.

Cite As

Jered Wells (2020). Simpson's 1/3 and 3/8 rules (https://www.mathworks.com/matlabcentral/fileexchange/33493-simpson-s-1-3-and-3-8-rules), MATLAB Central File Exchange. Retrieved .

Comments and Ratings (4)

Ariorh

Thanks for putting this up, it works great.

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.

MATLAB Release Compatibility
Created with R2008b
Compatible with any release
Platform Compatibility
Windows macOS Linux
Acknowledgements

Inspired by: Simpson's Rule Integration

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!