计算方法第二次上机作业.docxVIP

计算方法第二次上机作业 Gegebao 摘要：程序基于MATLAB，包括问题陈述、算法与程序、结果与分析、讨论四个部分。 一 、问题陈述 数学上已经证明了：0141+ 分别使用矩形、梯形和Simpson复合求积公式计算π的近似值。选择不同的h，对于每种求积公式，试将误差刻画成h的函数，并比较各方面的精度。是否存在某个h值，当低于这个值之后，再继续减少h的值，计算不再有所改进？为什么？ 实现Romberg求积方法，并重复上面的计算。 使用自适应求积方法重复上面的计算。 算法与程序 1．矩形求积法 考虑到对于具体的某一个h，不一定能整数个地覆盖积分空间，将会严重地影响算法精度。于是我们用n代替h，即将积分区间划分为n个区间，h=1/n。函数将输出积分结果Rec和这个结果对于精确的π的误差的对数值r。 function [Rec r ] = rec( n ) h=1/n; %求出对应的h Rec=0; for i=1:n x0=h*(i-1/2); %求出每一区间的中心点的横坐标 Rec=Rec+h*4/(1+x0^2); end r=log10(abs(Rec-pi)); %取其误差的对数 end 2.梯形求积法 和矩形法一样，采用n来标度取点密度。同样输出积分结果Lad和这个结果对于精确的π的误差的对数值r。 function [Lad r ] = lad( n) h=1/n; Lad=0; for i=1:n x0=h*(i-1); x1=x0+h; 取得每个区间的端点的横坐标 Lad=Lad+h/2*(4/(1+x0^2)+4/(1+x1^2)); end r=log10(abs(Lad-pi)); end 3.Simpson求积法 输入取点数目n，h=1/n, 输出积分结果Simp和这个结果对于精确的π的误差的对数值r。 function [Simp ,r ] = simp( n) h=1/n; Simp=0; for i=1:n x0=h*(i-1); x1=x0+h/2; x2=x1+h/2; Simp=Simp+h/6*(4/(1+x0^2)+4*4/(1+x1^2)+4/(1+x2^2)); end r=log10(abs(Simp-pi)); end 4.Romberg求积法 输入最初取点数目n（h=1/n）和需求递推下去的步数k, 输出积分结果Rom和这个结果对于精确的π的误差的对数值r。 function [Rom,r] = rom(n,k) for i=1:k; [a(i,1) c]=lad(2^(i-1)*n); %调用lad函数，求出递推矩阵第一列的数值 end for i=2:k for j=1:k-i+1 a(j,i)=(a(j+1,i-1)-4^(1-i)*a(j,i-1))/(1-4^(1-i)); end end Rom=a(1 ,k); r=log10(abs(Rom-pi)); end 5.分别基于前几种算法的自适应算法 （1）基于矩形法 输入初始的取点数n（h=1/n）和所需要的精度∈的对数值e，输出积分结果对于精确的π的误差的对数值r。 function [r] = Arec(n,e) temp0=0;temp1=1; %temp0为T（h），Temp1为T（h/2） while(log10(abs(temp0-temp1))e) temp0=temp1; n=n*2; temp1=rec(n); end r=log10(abs(temp1-pi)); end （2）基于梯形法 输入初始的取点数n（h=1/n）和所需要的精度∈的对数值e function [r] = Alad( n,e) temp0=0;temp1=1; while(log10(abs(temp0-temp1))e) temp0=temp1; n=n*2; temp1=lad(n); end r=log10(abs(temp1-pi)); end （3）基于Simpson法 输入初始的取点数n（h=1/n）和所需要的精度∈的对数值e function [r] = Asimp( n,e) temp0=0;temp1=1; while(log10(abs(temp0-temp1))e) temp0=temp1; n=n*2; temp1=simp(n); en