不知道要选什么的被积函数比较合适,
format long
fun=@(x)cos(x);
x0=2;
x1=3;
tol=10.^[-3:-2:-7,-8];
n=length(tol);
info=num2str(tol');
info=mat2cell(info,ones(n,1),size(info,2));
y=sin(x1)-sin(x0); %精确的解
for i=1:n
subplot(2,2,i);
[y1,n1,yiter1]=RectangInt(fun,x0,x1,tol);
[y2,n2,yiter2]=TrapzInt(fun,x0,x1,tol);
[y3,n3,yiter3]=SimpsonInt(fun,x0,x1,tol);
[y4,n4,yiter4]=RombergInt(fun,x0,x1,tol);
plot(1:log2(n1)+1,yiter1,1:log2(n2)+1,yiter2,1:log2(n3),yiter3,...
1:log2(n4)+1,yiter4([(1:log2(n4)+1).*(2:log2(n4)+2)/2]));
legend('矩形法','梯形法','Simpson','Romberg')
xlabel('对半次数')
ylabel('int');
title(info{i})
end
%%%%指定容限比较
abs([y1,y2,y3,y4]-y) %残差
[n1,n2,n3,n4]
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%矩形法是比较差的一种方法
fun=@(x)exp(-x.^2);
x0=0.5;
x1=2;
tol=1e-12;
y=erf(x1)-erf(x0);
tic; [y2,n2,yiter2]=TrapzInt(fun,x0,x1,tol); t2=toc;
tic; [y3,n3,yiter3]=SimpsonInt(fun,x0,x1,tol); t3=toc;
tic; [y4,n4,yiter4]=RombergInt(fun,x0,x1,tol); t4=toc;
figure
plot(1:log2(n2)+1,yiter2,1:log2(n3),yiter3,...
1:log2(n4)+1,yiter4([(1:log2(n4)+1).*(2:log2(n4)+2)/2]));
abs([y2,y3,y4]-y) %残差
[n2,n3,n4]
[t2,t3,t4]
legend('梯形法','Simpson','Romberg')
format
复制代码
在本例中。
对于1e-9容限,矩形法就比较吃力。
对于1e-15容限,梯形法也比较吃力。
|