%下面除了龙贝格外,其他均以此fun函数作为被积函数
%梯形、辛普森、科特斯:都是已知积分上下限和分段数,求解积分近似值。
function a = fun(x)
a = cos(x*x);
end
一、复化梯形:
function Tn = ComplexTrap(xL, xR, n)
% xL,xR:积分区间左右端点,n:分段数
% Tn :复化梯形积分结果
h = (xR-xL)/n;%求步长
Tn = 0;
x = xL:h:xR;
for k = 1:n
Tn = Tn + fun(x(k)) + fun(x(k+1));
end
Tn = Tn*h/2;
Tn = vpa(Tn,6);
end
二、复化辛普森:
function Sn = ComplexSimpson(xL, xR, n)
% xL,xR:积分区间左右端点,n:分段数
% Sn :复化辛普森积分结果
h = (xR-xL)/n;%求步长
Sn = 0;
x = xL:h:xR;
for k = 2:n+1
Sn = Sn + 4*fun((x(k)-0.5)) + 2*fun(x(k));
end
Sn = Sn + fun(xL) - fun(xR);
Sn = vpa(Sn,7);
Sn = Sn*h/6;
Sn = vpa(Sn,6);
end
三、复化科斯特:
function Cn = ComplexCotes(xL, xR,n)
% 判断积分区间个数是否是2的倍数,满足则进行计算,否则打印提示
if mod(n,4) ~= 0
print('等分区间数错误!');return
else
% 获取步长h
h = (xR-xL)/n;
%累积计算
Cn = 0;
for i = 1:n-3
Cn = Cn + 7*fun(xL+h*(i-1))...
+32*fun(xL+h*i)...
+12*fun(xL+h*(i+1))...
+32*fun(xL+h*(i+2))...
+7*fun(xL+h*(i+3));
end % 循环结束
Cn = Cn*h/90;
Cn = vpa(Cn,6);
end % if判断结束
end % 函数结束
四、龙贝格:
function I = romberg(fun,a,b,e)
% 使用龙贝格(Romberg数值求解公式)
% 输入例如:I=romberg(@(x)x^(3/2),0,1,0.000001)
% 判断输入参数是否足够
if nargin~=4
error('请输入需要求积分的f、上界和下界以及误差e')
end
k=0; % 迭代次数
n=1; % 区间划分个数
h=b-a; %上下界间距
T(1,1)=h/2*(fun(a)+fun(b));
d=b-a; %误差间距
while e<=d
k=k+1;
h=h/2;
sum=0;
% 计算第一列T
for i=1:n
sum=sum+fun(a+(2*i-1)*h);
end
T(k+1,1)=T(k)/2+h*sum;
% 迭代
for j=1:k
T(k+1,j+1)=T(k+1,j)+(T(k+1,j)-T(k,j))/(4^j-1);
end
n=n*2;
d=abs(T(k+1,k+1)-T(k,k));
end
T
I=T(k+1,k+1);