-
Notifications
You must be signed in to change notification settings - Fork 0
/
NumIntgrl.m
50 lines (45 loc) · 1.63 KB
/
NumIntgrl.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
function [method, result] = NumIntgrl(x, y)
% Numerical Integration utilizing Simpson's Rule. Defaults to the
% Trapezoidal Rule if Simpson's Rule cannot be applied.
%
% https://web.engr.oregonstate.edu/~webbky/MAE4020_5020_files/Section%208%20Integration.pdf
narginchk(2,2) % check # of input arguments
Nseg = length(y) - 1;
dx = ( max(x) - min(x) )/( length(x) - 1 ); % x-spacing
if mod(Nseg, 2) == 0 % Simpson's 1/3 rule: number of segments must be even
method = 'Simp 1/3';
sumOdd = 0; % sum over odd segments
for i=1:2:Nseg-1
sumOdd = sumOdd + y(i+1);
end
sumEven = 0; % sum over even segments
for i=2:2:Nseg-2
sumEven = sumEven + y(i+1);
end
result = 1/3*dx*(y(1) + 4*sumOdd + 2*sumEven + y(Nseg+1));
elseif mod(Nseg, 3) == 0 % Simpson's 3/8 rule: number of segments must be divisible by 3
method = 'Simp 3/8';
sum1 = 0;
for i=1:3:Nseg-2
sum1 = sum1+ y(i+1);
end
sum2= 0;
for i=2:3:Nseg-1
sum2 = sum2 + y(i+1);
end
sum3= 0;
for i=3:3:Nseg-3
sum3 = sum3 + y(i+1);
end
result = 3/8*dx*(y(1) + 3*sum1 + 3*sum2 + 2*sum3 + y(Nseg+1));
else % default to trapezoidal rule
method = 'Trap';
sum = 0;
for i=1:Nseg+1
sum = sum + y(i);
end
correction = ( y(1) + y(Nseg+1) )/2;
result = dx*(sum - correction);
% result = trapz(x, y); % compare against built-in MATLAB function
end
end