计算方法第二次上机作业

1、计算方法第二次上机作业Gegebao摘要：程序基于MATLAB，包括问题陈述、算法与程序、结果与分析、讨论四个部分。一 、问题陈述 数学上已经证明了：0141+x2dx=成立，所以可以通过积分来计算的近似值。（1） 分别使用矩形、梯形和Simpson复合求积公式计算的近似值。选择不同的h，对于每种求积公式，试将误差刻画成h的函数，并比较各方面的精度。是否存在某个h值，当低于这个值之后，再继续减少h的值，计算不再有所改进？为什么？（2） 实现Romberg求积方法，并重复上面的计算。（3） 使用自适应求积方法重复上面的计算。二、 算法与程序1矩形求积法考虑到对于具体的某一个h，不一定能整数个地覆

2、盖积分空间，将会严重地影响算法精度。于是我们用n代替h，即将积分区间划分为n个区间，h=1/n。函数将输出积分结果Rec和这个结果对于精确的的误差的对数值r。function Rec r = rec( n ) h=1/n; %求出对应的hRec=0;for i=1:n x0=h*(i-1/2); %求出每一区间的中心点的横坐标 Rec=Rec+h*4/(1+x02);endr=log10(abs(Rec-pi); %取其误差的对数end2.梯形求积法 和矩形法一样，采用n来标度取点密度。同样输出积分结果Lad和这个结果对于精确的的误差的对数值r。function Lad r = lad( n)

3、 h=1/n;Lad=0;for i=1:n x0=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:n x0=h*(i-1); x1=x0+h/2; x2=x1+h/2; Simp=Simp+h/6*(4/(1+x02)+4*4/(1+x12)+4/(1

4、+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: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); endendRom=a(1 ,k);r=log10(abs(Rom-pi);end5.

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);endr=log10(abs(temp1-pi);end （2）基于梯形法 输入初始的取点数n（h=1/n）和所需要的精度的对数值efunction r = Alad( n,e)temp0=

6、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:200 h0(i)=-i*7/200; n=fix(10(-h0

8、(i); r Rrec(i)=rec(n); endplot(h0,Rrec,'b');分析：求积方法的最高精度可以达到10(-14)到10(-15)之间，而当h取到10(-6)时，继续提高取点密度对算法精度没有改进。而在此之前，精度和h的关系几乎是线性的。2.梯形求积法运行代码：for i=1:200 h0(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取小开始波动，继续

9、减小h并不能保证精度的增加。而在此之前，精度和h的关系几乎是线性的。3.Simpson求积法运行代码：for i=1:200 h0(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所给出的值，计算机得到的误差为零，求对数运算无意义，故出现间

10、断点。4.三种算法的精度比较在进行前三个积分法的运算的时候，已经分别记录了结果，直接输出就行。运行代码：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:200 h0(

11、i)=-i*7/200; n=fix(10(-h0(i) r Rrom(i)=rom(n,k);endplot(h0,Rrom,'b');分析：可以看出，积分精度随h减小提高地非常快。在（-1,0）这个区间里，便出现了间断点，说明算法的精度已经达到MATLAB所能达到的精度，从而误差变为零，求对数运算无意义，出现间断点。纵观函数，最高精度应在10(-15)到10(-16)之间。6.自适应算法（取自适应精度为10(-10)）：（1）基于矩形算法 运行代码：for i=1:200 h0(i)=-i*7/200; n=fix(10(-h0(i) Rarec(i)=Arec(n,-10

12、);endplot(h0,Rarec,'b');(2)基于梯形算法：运行代码：for i=1:200 h0(i)=-i*7/200; n=fix(10(-h0(i) Ralad(i)=Alad(n,-10);endplot(h0,Ralad,'b');（3）基于Simpson算法 运行代码：for i=1:200 h0(i)=-i*7/200; n=fix(10(-h0(i) Rasimp(i)=Asimp(n,-10);endplot(h0,Rasimp,'b');（4）基于Romberg算法（ 固定递推步数为五步） 运行代码：for i=1:

13、200 h0(i)=-i*7/200; n=fix(10(-h0(i) Rarom(i)=Arom(n,-10);endplot(h0,Rarom,'b');分析：当本来的算法精度低于所给出的自适应精度时，自适应算法对本来算法的精度是有改进作用的，并使算法精度略高于所给出的自适应精度。但是当算法对于具体的h得到的结果的精度高于自适应给出的精度时，自适应算法则对原本的算法无改进作用。而且有些算法的精度随h的减小提高得特别快（比如说Romberg算法），使用自适应算法后算法精度会明显高于给出的自适应精度。四、讨论 算法的精度主要取决于算法本身、h的大小和机器误差。几种算法的精确程度不同，但都要受制于机器本身的误差，不可能无限地提高精度。当然，并不是说