资源描述:
《精品[精品]00数值分析课程设计》由会员上传分享,免费在线阅读,更多相关内容在应用文档-天天文库。
1、数值分析课程设计报告学号10801010223姓名王艳荣学号10801010208姓名何秋健学号10801010108姓名金绍荣学号10801010123姓名王松2011年6月2312数学与统计学院专业课程设计成绩评定书设计题目:Newton-Cotes和Romberg求积分信息与计算科学专业指导教师刘瑞华指导教师评语成绩:指导教师时间:答辩小组意见设计成绩:答辩组长:审定系主任:12题目:用熟悉的计算机语言编程上机完成(1)用Newton-Cotes公式计算积分的近似值,自己设置不同精度要求,对结
2、果进行比较分析。(2)用Romberg积分法计算积分的近似值,自己设置不同精度要求,对结果进行比较分析;与(1)的结果进行比较分析,谈谈你的体会。(3)记,在上面的计算中只取4位有效数字或7位有效数字,计算结果有什么不同。(4)上面计算精度可达8-20位有效数字吗?若可以请说明实现过程,并举例。一.摘要:在matlab环境下熟悉的运用计算机编程语言并结合Newton-Cotes和Romberg的理论基础对函数求积分,在运行完程序后以及对运行结果做出各方面的分析和比较。二.实验设计目地用熟悉的计算机语
3、言编程上机完成Newton-Cotes公式计算积分和Romberg积分法计算积分。三.Newton-Cotes和Romberg的理论基础牛顿柯斯特公式:设将积分区间[a,b]划分为n等份,步长h=(b-a)/n,选取等距节点=a+kh构造出的插值型求积分公式I,称为牛顿—柯斯特(Newton-Cotes)公式,式中称为柯特斯系数,引进变换x=a+th,则有龙贝格求积公式:(1).梯形法的递推化12设将区间[a,b]分为n等份,共有n+1个分点,如果将求积区间再二分一次,则分点增至2n+1个,我们将二
4、分前后两人积分值联系起来加以考察,注意到每个子区间[]经过二分只增加了一个分点,用复合梯形公式求得该子区间上的积分值为,注意,这里h=(b-a)/n代表二分前的步长,将每个子区间上的积分值相加得,从而利用上式可导出下列递推分式:(2).外推技巧从梯形公式出发,将区间[a,b]逐次二分可提高求积公式精度,当[a,b]分为n等份I若记,当区间[a,b]分为2n等份时,则有,并且有I.(3).龙贝格算法将上述外推技巧得到的公式重新引入记号等从而可将上述公式定成统一形式经过m(m=1.2.3……)次加速后,
5、余项便取下列形式:这方法通常称为理查森外推加速方法。设以表示二分次后求得的梯形值,且以表示序列{}的m次加速值,则依递推公式可得公式也称为龙贝格算法,计算过程如下:ⅰ.取k=0,h=b-a,求令(k记区间[a,b]二分次数)12ⅱ.求梯形值,即按递推公式计算ⅲ.求加速值,按公式a逐个求出T表的第k行其余各元素T(j=1,2,…,k).iv.若
6、T—T
7、<(预先设定的精度),则终止计算,并取TI;否则令k+1k转(2)继续计算.kh…0b-a…1①…2②③…3④⑤⑥…4⑦⑧⑨⑩……T表表指出了计算过程
8、,第2列h=给出了子区间长度,i表示第i步外推.可以证明,如果f(x)充分光滑,那么T表每一列的元素及对角线元素均收敛到所要求的积分值I.四.程序代码及运行结果牛顿柯斯特求积分:functiony=mulNewtonCotes(a,b,m,n)fun=@(x)exp(-x.^2);12xk=linspace(a,b,m+1)fori=1:ms(i)=NewtonCotes(fun,xk(i),xk(i+1),n);endy=sum(s);function[y,Ck,Ak]=NewtonCotes(f
9、un,a,b,n)%fun=@(x)exp(-x.^2)ifnargin==1[mm,nn]=size(fun);ifmm>=8errorendxk=fun(1,:);fk=fun(2,:);a=min(xk);b=max(xk);n=mm-1;endifnargin==4xk=linspace(a,b,n+1);ifisa(fun,'function_handle')fx=fun(xk);endendCk=cotescoeff(n);Ak=(b-a)*Ck;y=Ak*fx';12functionC
10、k=cotescoeff(n)fori=1:n+1k=i-1;Ck(i)=(-1)^(n-k)/factorial(k)/factorial(n-k)/n*quadl(@(t)intfun(t,n,k),0,n);endfunctionf=intfun(t,n,k)f=1;fori=[0:k-1,k+1:n]f=f.*(t-i);end运行结果>>y=mulNewtonCotes(0,1,2,1)xk=00.50001.0000y=0.7314>>vpa(y,7)an