function T=simpson(f_name, a, b, n)
% f_name为要求的定函数y=f(x)所在的程序文件名
% a为积分下限
% b为积分上限
% n为积分区间[a,b]划分成小区间的等份数
h = (b-a) / n;
x = a + (0:n) * h;
f = feval(f_name, x);
N = length(f) - 1;
if N == 1
fprintf('Data has only one interval')
return;
end
if N == 2
T = h / 3 * (f(1) + 4 * f(2) + f(3));
return;
end
if N == 3
T = 3 / 8 * h * (f(1) + 3 * f(2) + 3 * f(3) + f(4));
return;
end
T = 0;
if 2 * floor(N / 2) == N
T = h / 3 * (2 * f(N-2) + 2 * f(N-1) + 4 * f(N) + f(N+1));
m = N - 3;
else
m = N;
end
T = T + (h/3) * (f(1) + 4 * sum(f(2:2:m)) + 2 * f(m+1));
if m > 2
T = T + (h/3) * 2 * sum(f(3:2:m));
end
function y=simpson_f(x)
y = sin(x)
end
T = simpson('simpson_f', 0, pi, 4)
%=> y = 0 0.7071 1.0000 0.7071 0.0000
%=> T = 2.0046
function T=simpson(f_name, a, b, n)% f_name为要求的定函数y=f(x)所在的程序文件名% a为积分下限% b为积分上限% n为积分区间[a,b]划分成小区间的等份数h = (b-a) / n;x = a + (0:n) * h;f = feval(f_name, x);N = length(f) - 1;if N == 1 fprintf('Data has only one interval') return;endif N