《数值分析》课程设计报告范文

课程设计报告

课程设计题目：

学生姓名: 专 业 班 级 指导教师

年 月 日

题目：

用熟悉的计算机语言编程上机完成

（1）用Newton-Cotes公式计算积分exdx的近似值，自己设置不同精度要求，

012对结果进行比较分析。

（2）用Romberg积分法计算积分exdx的近似值，自己设置不同精度要求，对

012结果进行比较分析；与（1）的结果进行比较分析，谈谈你的体会。

（3）记f(x)ex，在上面的计算中f(x)只取4位有效数字或7位有效数字，计算结果有什么不同。

（4）上面计算精度可达8-20位有效数字吗？若可以请说明实现过程，并举例。

2一、摘要

在matlab环境下熟悉的运用计算机编程语言并结合Newton-Cotes和Romberg的理论基础对函数求积分，在运行完程序后以及对运行结果做出各方面的分析和比较。

二、设计目的

用熟悉的计算机语言编程上机完成Newton-Cotes公式计算积分和

Romberg积分法计算积分。

二、理论基础

牛顿柯斯特公式：设将积分区间[a,b]划分为n等份，步长h=（b-a）/n,

nn选取等距节点xk=a+kh构造出的插值型求积分公式I(baCk)f(k)，称为

k0(n)牛顿—柯斯特（Newton-Cotes）公式，式中Ck称为柯特斯系数，引进变换x=a+th,则有Ckh/(ba)龙贝格求积公式：

（1）.梯形法的递推化 设将区间[a,b]分为n等份，共有n+1个分点，如果将求积区间再二分一次，则分点增至2n+1个，我们将二分前后两人积分值联系起来加以考察，注意到每个子区间[

(n)n0n(n)j0jk(tj)(kj)dtnk!(nk)!0(1)nknn(tj)dt

j0jkx,xkk1]经过二分只增加了一个分点

xh4k1212(xkxk1)，用复合梯形公式求得该子区间上的积分值为

=(b-a)/n代表二分前的步长，将每[f(x)2f(xk1)f(xk1)]，注意，这里h

2个子区间上的积分值相加得T2nh4n1[f(xk0k)f(xk1)]n1h2n1k0f(xk1)2，从而利

用上式可导出下列递推分式：T12nT2nh2k0f(xk1)

2（2）.外推技巧 从梯形公式出发，将区间[a,b]逐次二分可提高求积公式精度，当[a,b]分为n等份ITnba12hf2n(),[a,b],hban

，并且有

若记TnT(h)，当区间[a,b]分为2n等份时，则有TT(h)Iba122nhT()2hf2''(),limT(h)T(0)I.

h0(3).龙贝格算法 将上述外推技巧得到的公式重新引入记号

TT0(h)T(h),(h)C(h),mTT1(h)S(h)

(h)R(h) 等从而可将上述公式定成统一形式

23T(h)m4T14mh()m1214mT1m1(h)经过m(m=1.2.3……)次加速后，余项便

取下列形式：

T(h)Im1h2(m1)2h2(m1)........

这方法通常称为理查森外推加速方法。 设以T(k)m表示二分次后求得的梯形值，且以T(k)m(k)m表示序列｛T,k1.2......(k)m｝的m次加速值，

则依递推公式可得T4T14mm(k1)m114mT1(k)m1

公式也称为龙贝格算法，计算过程如下： ⅰ.取k=0,h=b-a,求T(0)0h2[f(a)f(b)]

令1k（k记区间[a,b]二分次数）

ⅱ.求梯形值T0(ba2k)，即按递推公式计算T(k)0

(kj)jⅲ.求加速值，按公式a逐个求出T表的第k行其余各元素Tiv.若|T

(0)k（j=1,2,…,k）.

(0)—T

(0)k1|<(预先设定的精度)，则终止计算，并取TkI；否则令

k+1k转（2）继续计算. k 0 1 2 3 4  h b-a ba2ba4ba8ba16T 0(k)T 1(k)T (k)2 T (k)3 T (k)4 … … … … … T 0(1)0(0) (0) T① T 1(2)0(1)1 (0)2T② T③ T(3)0(2)1 (0)3T④ T⑤ T⑥ T2(1) (0)4T⑦ T⑧ T01(4)(3)(2)2⑨ T⑩ T3(1) … …       T 表 bak表指出了计算过程，第2列h=2给出了子区间长度，i表示第i步外推.可以

证明，如果f(x)充分光滑，那么T表每一列的元素及对角线元素均收敛到所要求的积分值I.

三、程序代码及运算结果

牛顿柯斯特求积分：

function y=mulNewtonCotes(a,b,m,n) fun=@(x)exp(-x.^2); xk=linspace(a,b,m+1) for i=1:m

s(i)=NewtonCotes(fun,xk(i),xk(i+1),n);

end y=sum(s);

function [y,Ck,Ak]=NewtonCotes(fun,a,b,n) %fun=@(x)exp(-x.^2) if nargin==1

[mm,nn]=size(fun); if mm>=8 error end xk=fun(1,:); fk=fun(2,:); a=min(xk); b=max(xk); n=mm-1; end

if nargin==4

xk=linspace(a,b,n+1);

if isa(fun,'function_handle') fx=fun(xk); end end

Ck=cotescoeff(n); Ak=(b-a)*Ck; y=Ak*fx';

function Ck=cotescoeff(n) for i=1:n+1 k=i-1;

Ck(i)=(-1)^(n-k)/factorial(k)/factorial(n-k)/n*quadl(@(t)intfun(t,n,k),0,n);

end

function f=intfun(t,n,k) f=1;

for i=[0:k-1,k+1:n] f=f.*(t-i);

end

运行结果

>> y=mulNewtonCotes(0,1,2,1) xk =

0 0.5000 1.0000 y = 0.7314 >> vpa(y,7) ans = .7313703

>> y=mulNewtonCotes(0,1,5,2) xk =

0 0.2000 0.4000 y = 0.7468 >> vpa(y,7) ans = .746824 >>

龙贝格求积分

function z=romberg(a,b,e) h=(b-a);

f=@(x)exp(-x.^2);

TT(1,1)=h.*(f(b)+f(a))/2;

0.6000 0.8000 1.0000 k=2;

TT(1,2)=TT(1,1)./2+h/2.*f(a+h/2); TT(2,1)=TT(1,2).*4/3-TT(1,1)./3; z=TT(2,1);

while abs((TT(k,1)-TT(k-1,1))./TT(k,1))>=e k=k+1;h=h./2; for j2=1:2.^(k-2)

ff(1,j2)=f(a+h*(j2-1/2)); end

fff=sum(ff).*h/2; TT(1,k)=TT(1,k-1)./2+fff for j1=2:k

TT(j1,k-j1+1)=4^(j1-1).*TT(j1-1,k-j1+2)./(4^(j1-1)-1)-TT(j1-1,k-j1+1)/(4^(j1-1)-1); z=TT(j1,k-j1+1); end end

function r=f() r=@(x)exp(-x.^2)

运算结果

>> z=romberg(0,1,10e-2) z = 0.7472 >> vpa(z,7) ans = .7471804

>> z=romberg(0,1,10e-5) TT =

0.6839 0.7314 0.7430 0.7472 0 0 TT =

0.6839 0.7314 0.7430 0.7459 0.7472 0.7469 0 0 0.7468 0 0 0 z = 0.7468 >> vpa(z,7) ans = .7468240

>> z=romberg(0,1,10e-10) TT =

0.6839 0.7314 0.7430 0.7472 0 0 TT =

0.6839 0.7314 0.7430 0.7472 0.7469 0 0.7468 0 0 TT =

0.6839 0.7314 0.7430 0.7472 0.7469 0.7468 0.7468 0.7468 0 0.7468 0 0 TT =

0.6839 0.7314 0.7430 0.7472 0.7469 0.7468 0.7468 0.7468 0.7468 0.7468 0.7468 0 0.7459 0 0 0.7459 0 0 0 0.7459 0.7468 0 0 0.7466 0 0 0 0.7466 0 0 0 0.7468 0 0 0

0.7468 0 0 0 0 0 z = 0.7468 >> vpa(z,7) ans = .7468241

四、结果分析

根据牛顿柯斯特求积分的原理将所求函数的积分区间分为多个小的积分

区间，先求出每个小积分区间上的函数值，然后再将每个小区间上的求积结果加起来就是我们所要求的总函数的积分值，当函数区间所分的小区间的个数越多的时候总函数所计算出来的结果就精确，其原因就是所分的区间数越多，计算时每个小区所带来的误差就越小，当将总的积分值加起来的时候所带来的总的误差也就越小，所以最后的结果的精度越高，而用龙贝格求积和牛顿柯斯特也一样是要将总区间分为很多小的相等的区间，只是他们所用的计算原理不一样，当用此方法求积分时，所设的误差限越小，所求得的积分值就越是精确。

五、设计心得 自己完成。 六、参考文献

[1]孙祥，徐流美，吴清。Matlabl 7.0

版社，2005.

[2]薛毅。数值分析与实验【M】.北京：北京理工大学出版社，2005. [3]汪卉琴，刘目楼。数值分析【M】。北京：冶金工业出版社，2004.

[4]丁丽娟，程杞元。数值计算方法【M】。北京：北京理工大学出版社，2005. [5]肖伟，刘忠，曾新勇，等。Matlabl程序设计与应用【M】。北京：清华大学出版社，2005.

[6]李庆扬，易大义，王能超。现代数值分析。北京：高等教育出版社，1995.

基础教程【M】。北京：清华大学出