RES = SIMPSON(Y) computes an approximation of the integral of Y via
辛普森的1/3俄文le (with unit spacing). Simpson's 1/3 rule uses
quadratic interpolants for numerical integration. To compute the
积分间距不同,相乘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
辛普森的1/3俄文le. 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.