version 1.2.0.0 (3.11 KB) by
Jered Wells

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

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.

Jered Wells (2021). 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 .

Created with
R2008b

Compatible with any release

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

Start Hunting!Create scripts with code, output, and formatted text in a single executable document.

AriorhThanks for putting this up, it works great.

Evan Provaneric rossAndrew DavisThanks 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.