1、计算方法第二次上机作业Gegebao摘要：程序基于 MATLAB，包括问题陈述、算法与程序、结果与分析、讨论四个部分。一 、问题陈述 数学上已经证明了： 成立，所以可以通过积分来计算 的近似10 41+2= 值。（1） 分别使用矩形、梯形和 Simpson 复合求积公式计算 的近似值。选择不同的 h，对于每种求积公式，试将误差刻画成 h 的函数，并比较各方面的精度。是否存在某个 h 值，当低于这个值之后，再继续减少 h 的值，计算不再有所改进？为什么？（2） 实现 Romberg 求积方法，并重复上面的计算。（3） 使用自适应求积方法重复上面的计算。二、 算法与程序1矩形求积法考虑到对于具体的
2、某一个 h，不一定能整数个地覆盖积分空间，将会严重地影响算法精度。于是我们用 n 代替 h，即将积分区间划分为 n 个区间，h=1/n。函数将输出积分结果 Rec 和这个结果对于精确的 的误差的对数值 r。function Rec r = rec( n ) h=1/n; %求出对应的hRec=0;for i=1:nx0=h*(i-1/2); %求出每一区间的中心点的横坐标Rec=Rec+h*4/(1+x02);endr=log10(abs(Rec-pi); %取其误差的对数end2.梯形求积法和矩形法一样，采用n来标度取点密度。同样输出积分结果Lad和这个结果对于精确的 的误差的对数值r。fu
3、nction Lad r = lad( n) h=1/n;Lad=0;for i=1:nx0=h*(i-1);x1=x0+h; 取得每个区间的端点的横坐标Lad=Lad+h/2*(4/(1+x02)+4/(1+x12);endr=log10(abs(Lad-pi);end3.Simpson求积法输入取点数目n ，h=1/n, 输出积分结果Simp 和这个结果对于精确的 的误差的对数值r。function Simp ,r = simp( n)h=1/n; Simp=0;for i=1:nx0=h*(i-1);x1=x0+h/2;x2=x1+h/2;Simp=Simp+h/6*(4/(1+x02)
4、+4*4/(1+x12)+4/(1+x22);endr=log10(abs(Simp-pi);end4.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函数，求出递推矩阵第一列的数值endfor i=2:kfor j=1:k-i+1a(j,i)=(a(j+1,i-1)-4(1-i)*a(j,i-1)/(1-4(1-i);endendRom=a(1 ,k);r=log10(ab
5、s(Rom-pi);end5.分别基于前几种算法的自适应算法（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);endr=log10(abs(temp1-pi);end（2）基于梯形法输入初始的取点数n（h=1/n）和所需要的精度 的对数值efunction r = Alad( n
6、,e)temp0=0;temp1=1;while(log10(abs(temp0-temp1)e)temp0=temp1;n=n*2; temp1=lad(n);endr=log10(abs(temp1-pi);end（3）基于Simpson法输入初始的取点数n（h=1/n）和所需要的精度 的对数值efunction r = Asimp( n,e)temp0=0;temp1=1;while(log10(abs(temp0-temp1)e)temp0=temp1;n=n*2; temp1=simp(n);endr=log10(abs(temp1-pi);end（4）基于Romberg法输入初始的
7、取点数n（h=1/n）和所需要的精度 的对数值e，取递推步数k=5function r = Arom( n,e)temp0=0;temp1=1;while(log10(abs(temp0-temp1)e)temp0=temp1;n=n*2; temp1=rom(n,5);endr=log10(abs(temp1-pi);end三、结果与分析使用之前的积分函数进行运算。考虑到h的变化范围比较大，将h 等比例地从1取到10(-7),取200个点。最后将h和误差按对数输出成图像。1.矩形积分法运行代码：for i=1:200h0(i)=-i*7/200;n=fix(10(-h0(i);r Rrec(
8、i)=rec(n); endplot(h0,Rrec,b);分析：求积方法的最高精度可以达到10(-14) 到10(-15) 之间，而当h 取到10(-6)时，继续提高取点密度对算法精度没有改进。而在此之前，精度和h的关系几乎是线性的。2.梯形求积法运行代码：for i=1:200h0(i)=-i*7/200;n=fix(10(-h0(i);r Rlad(i)=lad(n);endplot(h0,Rlad,b);分析：求积方法的最高精度可以达到10(-14) 到10(-15) 之间，而当h 取到10(-6)时，精度随 h取小开始波动，继续减小 h并不能保证精度的增加。而在此之前，精度和h 的关
9、系几乎是线性的。3.Simpson求积法运行代码：for i=1:200h0(i)=-i*7/200;n=fix(10(-h0(i)r Rsimp(i)=simp(n);endplot(h0,Rsimp,b);分析：求积方法的最高精度可以达到10(-15) 到10(-16) 之间，而当h 取到10(-2)时，精度随 h取小开始波动，继续减小 h并不能保证精度的增加。而在此之前，精度和h 的关系几乎是线性的。图像上存在间断点，实际上是由于求积方法的精度已经超过MATLAB所给出的 值，计算机得到的误差为零，求对数运算无意义，故出现间断点。4.三种算法的精度比较在进行前三个积分法的运算的时候，已经
10、分别记录了结果，直接输出就行。运行代码：plot(h0,Rrec,b);hold on;plot(h0,Rlad,r);hold on;plot(h0,Rsimp,k);hold on;分析：对于相同的h 值，在很长的范围内， Simpson的精度明显高于另两种算法，而矩形法又比梯形法略优。三种算法的极限精度以Simpson为最好，可以达到10（-15）到10（-16）之间。5Romberg积分法（固定递推步数为 5，基于梯形算法）运行代码：k=5for i=1:200h0(i)=-i*7/200;n=fix(10(-h0(i)r Rrom(i)=rom(n,k);endplot(h0,Rro
11、m,b);分析：可以看出，积分精度随h减小提高地非常快。在（-1,0 ）这个区间里，便出现了间断点，说明算法的精度已经达到MATLAB所能达到的精度，从而误差变为零，求对数运算无意义，出现间断点。纵观函数，最高精度应在10(-15)到10(-16)之间。6.自适应算法（取自适应精度为 10(-10)）：（1）基于矩形算法运行代码：for i=1:200h0(i)=-i*7/200;n=fix(10(-h0(i)Rarec(i)=Arec(n,-10);endplot(h0,Rarec,b);(2)基于梯形算法：运行代码：for i=1:200h0(i)=-i*7/200;n=fix(10(-h
12、0(i)Ralad(i)=Alad(n,-10);endplot(h0,Ralad,b);（3）基于Simpson算法运行代码：for i=1:200h0(i)=-i*7/200;n=fix(10(-h0(i)Rasimp(i)=Asimp(n,-10);endplot(h0,Rasimp,b);（4）基于Romberg算法（ 固定递推步数为五步）运行代码：for i=1:200h0(i)=-i*7/200;n=fix(10(-h0(i)Rarom(i)=Arom(n,-10);endplot(h0,Rarom,b);分析：当本来的算法精度低于所给出的自适应精度时，自适应算法对本来算法的精度是
13、有改进作用的，并使算法精度略高于所给出的自适应精度。但是当算法对于具体的h 得到的结果的精度高于自适应给出的精度时，自适应算法则对原本的算法无改进作用。而且有些算法的精度随h的减小提高得特别快（比如说Romberg算法），使用自适应算法后算法精度会明显高于给出的自适应精度。四、讨论算法的精度主要取决于算法本身、h的大小和机器误差。几种算法的精确程度不同，但都要受制于机器本身的误差，不可能无限地提高精度。当然，并不是说算法本身的改进对于相同的机器精度没有意义。达到机器所允许的极限精度的时间空间代价对于不同算法是不一样的，寻求用更少的空间和时间达到更高的精度便是改进算法的意义所在。从另一个角度讲，机器精度本来就应该纳入算法的设计中。在运算过程中，减少机器误差的影响，从而使最终结果达到更高的精度，也是很重要的。不同算法在这方面还是有差异的。以本次作业为例，在相同的机器精度下，Simpson算法的极限精度便高于矩形算法和梯形算法。
