搜档网
当前位置:搜档网 › 计算方法第三次上机实习报告

计算方法第三次上机实习报告

计算方法第三次上机实习报告
计算方法第三次上机实习报告

实验报告

课程名称: 计算方法 指导老师: 太英 成绩:

实验名称: 第三次上机作业 实验类型: matlab 同组学生: 一、实验目的和要求(必填) 二、实验容和原理(必填) 三、主要仪器设备(必填) 四、操作方法和实验步骤 五、实验数据记录和处理

六、实验结果与分析(必填)

七、讨论、心得

一、实验目的

用龙贝格算法计算积分I =∫ssss s ss 01

,要求误差不超过ε=12

×105

二、实验原理

龙贝格算法是由递推算法得来的。由梯形公式得出辛普森公式得出柯特斯公式最后得到龙贝格公式。设将求积区间[a ,b]分为n 个等分,则一共有n+1个等分点,k x a kh =+,0,1,b a

h k n

-=

=,n 。这里用n T 表示复化梯形法求得的积分

值,其下标n 表示等分数。

先考察下一个字段[1,k k x x +],其中点()112

1

2

k k k x x x ++=

+,在该子段上二分前后两个积分值 ()()112

k k h

T f x f x +=

+???? ()()21124k k k h T f x f x f x ++?

???

=+

+?? ???????

显然有下列关系 2112122k h T T f x +?

?=+ ???

将关系式关于k 从0到n-1累加求和,即可得递推公式1

210

2122n n n k k h T T f

x -+=??

=+ ???

需要强调指出的是,上式中的b a h n -=

代表二分前的步长,而12

12k x a k h +?

?=++ ???

根据梯形法的误差公式,积分值n T 的截断误差大致与2h 成正比,因此步长减半后误差将减至四分之一,即有

211

14

n n T T -≈- 将上式移项整理,知 221

1()3

n n n T T T -≈-

按上式,积分值2n T 的误差大致等于21

()3

n n T T -,如果用这个误差值作为2n T 的

一种补偿,可以期望,所得的()222141

333

n n n n n T T T T T T =+

-=-应当是更好的结果。 组合得到的近似值T 直接验证,用梯形二分前后的两个积分值n T 和2n T 按式组合,

结果得到辛普森法的积分值241

33

n n n S T T =-

再考察辛普森法。其截断误差与4h 成正比。因此,若将步长折半,则误差相应的减至十六分之一。既有

21

16n n I S I S -≈- 由此得 21611515

n n I S S ≈

- 不难验证,上式右端的值其实就等于n C ,就是说,用辛普森法二分前后的两个积分值n S 和2n S ,在按上式再做线性组合,结果得到柯特斯法的积分值n C ,既有

2161

1515n n n C S S ≈

- 重复同样的手续,依据斯科特法的误差公式可进一步导出龙贝格公式

2641

6363

n n n R C C ≈

- 在步长二分的过程中运用公式加工三次,就能将粗糙的积分值n T 逐步加工成精度较高的龙贝格n R ,或者说,将收敛缓慢的梯形值序列n T 加工成熟练迅速的龙贝

值序列n R ,这种加速方法称龙贝格算法。

三、实验过程

1.流程图

2.程序代码

clc

clear all;

format long

a=input('请输入你要求得积分的下限:');

b=input('请输入你要求得积分的上限:');

e=input('请输入你要求得积分的结束精度:');

k=input('请输入你要求得积分的最大次数:');

fx=(x)sin(x)/x;

lbg(f,a,b,k,e)

function lbg(fx,a,b,k,e)

T=zeros(k,k);

T(1,1)=(b-a)*(1+fx(b))/2;

for i=1:k

m=0;

for j=1:2^(i-1)

m=m+fx(a+(2*j-1)*(b-a)/(2^i));

end

T(i+1,1)=0.5*T(i,1)+(b-a)*m/2^i;

for n=1:i

T(i+1,n+1)=(4^n*T(i+1,n)-T(i,n))/(4^n-1); end

if abs(T(i+1,i+1)-T(i,i))<=e & i>=4

break;

else

;

end

end

for i=1:k

if T(i,1)==0

j=i;

break;

else

;

end

end

if j==k

error('所求次数不够或不可积') else

;

end

T=T(1:j-1,1:j-1)

disp('所求的积分值为:')

disp(T(j-1,1))

end

function fx = f(x)

if x == 0

fx = 1;

else

fx = sin(x) / x;

end

end

3.实验结果

四、实验心得

这是计算方法的最后一次上机实验,经过了前几次在MATLAB的摸爬滚打之后,我对MATLAB的使用比一开始熟练了很多。这次实验之中对基本函数的使用和对函数的调用等等我甚至已经可以做到不用借助网络搜索寻找函数。(第一次实验基本上就是需要一个函数就上网搜索他的意思和使用方法)。通过对龙贝格方法的使用,我体会到了龙贝格算法的快速收敛性(计算基本是瞬间完成了)。这也正是龙贝格算法的最大优势。

相关主题