数值计算方法

数值计算方法比较同等计算精度要求下变步长梯形法、变步长Simpson法、Romberg积分法在计算定积分时的速度及原因分析指导老师：王继红学生：陈奇 学号：080712110057 专业：电子信息科学与技术 时间：2011 年11 月22 日摘要本文由一道习题出发，主要论述了在同等计算精度要求下变步长 梯形法、变步长Simpson法、Romberg积分法在计算定积分时的速度 差异，及其原因分析。并指出在不同条件下如何选择数值计算方法以 节约计算时间。关键词：变步长、梯形法、Simpson、Romberg、积分法。目录摘 要1第一章 引言31.1 概述31.2 求积公式31.2.1 机械求积公式31.2.2 插值型求积公式41.2.3 复合求积公式51.2.4 Romberg 积分法61.2.5 Gauss 求积公式6第二章 程序及运行结果72.1 题目：72.2 解法及过程72.2.1 变步长梯形法：72.2.2 变步长 Simpson 法数值积分92.2.3 Romberg 积分法 132.3 结果汇总 17第三章 结果分析 183.1 结果简析 183.2 原因探究 19第四章 结论 22参考文献 23第一章 引言1.1 概述求积分r f( x)dx是数学科学的中心课题之一，主要途径有两条：a1. 微积分数学基本定理主要运用 Newton-Leibniz 公式求积分，其关键是找原函数。但大量被积函数找不到用 初等函数表示的原函数，更不用说解析式更加复杂的函数和用函数表给出的函数了。Newton-Leibniz 公式能力有限！2. 数值积分数值积分的起源可以追溯到古代。曲直转化极限趋近或“割圆术”正是数值积分的基本 思想。目前已经提出了许多数值积分方法，随着计算机技术的迅猛发展，数值积分这类看起 来原始、粗糙的方法日愈受到人们的重视。目前， MATLAB 等著名数学软件都提供有数值 积分功能。计算机时代主要用数值积分公式求积分。1.2 求积公式1.2.1 机械求积公式1.中矩形公式、梯形公式、Simpson公式中矩形公式：Jb f (x) dxa梯形公式：J b f ( x ) dxaSimpson 公式(抛物线公式)2. 机械求积公式J b f ( x ) dxab 一 aa + bq (f (a) + 4 f () + f (b)62Jb f (x) dx Q 为 A f (x )kk ak=0称 x 为求积节点，kA 为求积系数，它们与被积函数 f (x) 无关3.代数精度定义1：若一个机械求积公式对f (x) = xj (j = 0,1,m)准确成立，但对f (x) = xm+1不 准确成立，就说它具有m次代数精度。1.2.2 插值型求积公式1.插值型求积公式f ( x )l ( x)dxkkJbf (x) dx a J bp( x) dx =f b工naaak 二 0厦J b l (x)dxf (x )kkak 二 0其中l (x)=kx 一 xj 二 0, j 丰 k kj定理1具有n+1个求积节点的机械求公式至少具有n次代数精度的充要条件是它是插 值型的。2. Newton-Cotes 公式A f (x ) = (b a)为 Cnf (x )。kkkkak 二 0k 二 03. Newton-Cotes 公式代数精度由定理1知，n阶Newton-Cotes公式至少具有n次代数精度。可以证明，n为偶数时n阶 Newton-Cotes 公式的代数精度至少为 n+1。4. Newton-Cotes 公式的余项为了刻画数值积分公式的精度，除了使用“代数精度”，常常还必须讨论余项(截断误差)。如不特别声明时：I = Jb f (x) dx，ab 一 a=(f (a) + f (b)，12b 一 aa + bS =(f (a) + 4 f () + f (b) o1 2 2称I 一 T为梯形公式的余项,1I - S为Simpson公式的余项。1定理 2 设f (x) G C2a,b，则(b 一 a)312f 忆),定理 3 设f (x) g C4a,b，则1 S1 =册心G W1.2.3 复合求积公式1.复合梯形公式 从几何上看，积分区间较大时，在整个区间上用梯形公式计算公式计算效果不好(见图1-1)，分成若干个小区间，在每个区间上用梯形公式计算，然后加起来，效果好(见图1-2) 这种分割求和的方法就是复合求积。从图形中可以看出分割越细，精度越高。a) + 2 Z f (x ) + f (b) = Tknk 二 1其中，h = 山为步长，称上式为复合梯形公式，称T为复合梯形积分值。nn2.复合 Simpson 公式Simpson 公式的数值积分公式为：Jb f (x) dx =aa + bf (a) + 4 f () + f (b)2它的几何意义为用通过三点(a, f (a),(出2, f (S), (b, f (b)的抛物线围成的曲22边形面积来代替给定的函数的积分，精度高一些的还有Simpson3.8公式:Jb f(x)dx =a)+ f (b)其中同梯形公式一样，也有复合 Simpson 公式:Jb f(x)dx =ah n 一1Z f (x ) + 4f (x _) + f (x )6kk + 1k +12k=02f(x )= f(a), f(x )= f(b) 0nx + x ” b 一 a x _ = k k+1 , h =k + 22n1.2.4 Romberg 积分法Romberg 积分法是用 Richardson 外推算法来加快复合梯形求积公式的收敛速度，它的算 法如下，其中T是一系列逼近原定积分的Romberg积分值。m(1)计算b 一 aT(。)= f (a) + f (b)m2(2)对k = 1,2,3，计算下列各步：1b 一 a 2k-iT(k) =Tk-1 + 工 f (a +1212k 一1j=1(2 j - 1)( b - a)2 k对 m = 1,2,k 和 i = k, k 一 1, k 一 2,，计算:T i 一1 =m +14 m T i T i -1mm 。4m 一 11.2.5 Gauss 求积公式考虑插值型求积公式J1 f (x) dx一1a 为 A f (x ),kkk = 1其中J一119 世纪初期 Gauss 证明了存在唯nx xjdxx 一 xj = 1, j 丰 k kj种选择求积节点的方法，使得插值型求积公式具有 2n-1 次代数精度，并且这是可能达到的最高代数精度。定义 2：当插值型求积公式具有 2n-1 次代数精度时，称其为 n 点 Gauss 求积公式，称其求积节点x , x ,x为Gauss点。1 2 n-1,1上的Gauss求积公式，那么怎样用这些公式求Jbf (x)呢？首先用变换 ab 一 a a + bx =1 +22 将a,b上的积分转化为-1,1上的积分bb 一 aJ f (x) dx =a21 b 一 a a + b J f ( t +) dt一12 2第二章 程序及运行结果2.1 题目：给定积分exdx和J31 dx ,分别用下列方法计算积分值，要求准确到10-5,并比较分 11 x析计算时间。1）变步长梯形法；2）变步长 Simpson 法；3）Romberg 积分法。2.2 解法及过程2.2.1 变步长梯形法：梯形法数值积分采用的梯形公式是最简单的数值积分公式，其计算表达式为：ba 一 bJ f (x) dx = f (a) + f (b)a2由于用梯形公式来求积分十分粗糙，误差也比较大，后来改进后提出了复合梯形公式：Jb f (x) dx = f (a) + 2 S f (x ) + f (b)a2kk 二 1其中”b 一 ah =。n在 MATLAB 中编程实现的复合梯形（变步长梯形）公式的函数为： CombineTraprl。 功能：复合梯形公式求函数的数值积分。调用格式：I,step=CombineTraprl（f,a,b,eps）。其中，f 为函数名；a 为积分下限；b 为积分上限；eps 为积分精度；I 为积分值；Step 为积分划分的区间个数2.2.1.1 程序function I,step,t = CombineTraprl(f,a,b,eps)%复合梯形公式求函数f在区间a,b上的定积分%函数名：f%积分下限： a% 积分上限： b% 积分精度： eps%积分值： I% 积分划分的区间个数： stept0=cputime;if(nargin=3)eps=1.0e-4;endn=1;h=(b-a)/2;I1=0;I2=(subs(sym(f),findsym(sym(f),a)+subs(sym(f),findsym(sym(f),b)/h; while abs(I2-I1)epsn=n+1;h=(b-a)/n;I1=I2;I2=0;for i=0:n-1x=a+h*i;x1=x+h;I2=I2+(h/2)*(subs(sym(f),findsym(sym(f),x)+.subs(sym(f),findsym(sym(f),x1);endendI=I2;step=n;t=cputime-t0;%计算程序运行时间2.2.1.2计算3 exdx积分值，精确到10-5：1在 MATLAB 命令窗口中输入下列命令： f,step,t=CombineTraprl(exp(x),1,3,1.0e-5)f =17.367770318178550step =106t =117.8431554000000单 MATLAB 7.10.0 (R201iFile Edit Debug Parallel Desktop Window 旦巳 Ip7*0 已 | 語芒)专 為晾韵 | Current Folder: D:UsersM SI$6 的文档MATLAE 十 口 宦)0 Whats NewShort cuts 回 How to ACurrent F仆 K XName *E mfft矗CombineTrapd.m 矗i IntSimpson.mi 吏j Lagarange.mi 划 Romberg.miCommand WindowCf-1 New to MATILAB? Watch this 甘id&p, see Derng or read GetxT-17.367770318178550st ep =106Workspace倒題聒 | | Select data t. Name Value田田田17.3678106117.8432Command HistoryDetailsMinMax17.367817.3678106106117.8432 117.8432vpaffj 7) j st epj t =Romb erg C 1/ (f j st epj t =Romberg C 1/ (x)J 3 lj 3 ：format long“tjffj st epj t=Romberg ( esp (x)1 “ fj stepj t=CombineTraprl exp (f Start图1在MATLAB中用变步长梯形法计算exdx积分值12.2.1.3计算J31 dx积分值，精确到10-5：1x在 MATLAB 命令窗口中输入下列命令： f,step,t=CombineTraprl(1/(x),1,3,1.0e-5)f =1.098797422474723step =40t =17.144509900000032.2.2 变步长 Simpson 法数值积分Jb f (x) dx =aSimpson 公式的数值积分公式为：a + bf (a) + 4 f () + f (b)2它的几何意义为用通过三点(a, f (a),(出2, f (S), (b, f (b)的抛物线围成的曲Jb f (x) dx =边形面积来代替给定的函数的积分，精度高一些的还有Simpson3f8公式:b 一 a2 a + ba + 2 bf (a) + 3 f () + 3 f () + f (b)图2在MATLAB中用变步长梯形法计算丄dx积分值 1x 同梯形公式一样，也有复合 Simpson 公式：pF (x) dx = S / (x ) + 4 / (x _) + / (x )a6KK + 1K +12K 二 0/ (x ) = / (A), / (x ) = / (B )0n其中x + x”B ax kK 1 , H K+22n在MATLAB中编程实现的Simpson系列公式的函数为：IntSimpson。 功能： Simpson 系列公式求函数的数值积分。调用格式： I,step=IntSimpson(f,a,b,type,eps) 。其中，f 为函数名；a 为积分下限；b 为积分上限；type 为 Simpson 公式的类型；eps 为积分精度；I 为积分值；step 为积分划分的区间个数。2.2.2.1 程序function I,step,t = IntSimpson(f,a,b,type,eps)%Simpson系列公式求函数f在区间a,b上的定积分%函数名：f%积分下限： a% 积分上限： b%Simpson公式的类型：type (可取1, 2, 3)% 积分精度： eps%积分划分的子区间个数： stept0=cputime;if(type=3 & nargin=4)disp(缺少参数！)；endI=0；switch typecase 1,%Simpson 公式I=(b-a)/6)*(subs(sym(f),findsym(sym(f),a)+. 4*subs(sym(f),findsym(sym(f),(a+b)/2)+. subs(sym(f),findsym(sym(f),b)；step=1；case 2,%Simpson3/8 公式I=(b-a)/8)*(subs(sym(f),findsym(sym(f),a)+. 3*subs(sym(f),findsym(sym(f),(2*a+b)/3)+. 3*subs(sym(f),findsym(sym(f),(a+2*b)/3)+. subs(sym(f),findsym(sym(f),b)；step=1；case 3,%复合 Simpson 公式n=2；h=(b-a)/2；I1=0；I2=(subs(sym(f),findsym(sym(f),a)+subs(sym(f),findsym(sym(f),b)/h； while abs(I2-I1)epsn=n+1；h=(b-a)/n；I1=I2；I2=0；for i=0:n-1x=a+h*i； x1=x+h；I2=I2+(h/6)*(subs(sym(f),findsym(f),x)+. 4*subs(sym(f),findsym(sym(f),(x+x1)/2)+. subs(sym(f),findsym(sym(f),x1)；endendI=I2；step=n；endt=cputime-t0;2222计算J3 exdx积分值，精确到10-5：1在 MATLAB 命令窗口中输入下列命令： f,step,t=IntSimpson(exp(x),1,3,3,1.0e-5) f =17.367269778942639step =9 t =1.528809800000005图3在MATLAB中用变步长Simpson法计算J3 exdx积分值12223计算f31 dx积分值，精确到10-5： 1x在 MATLAB 命令窗口中输入下列命令： f,step,t=IntSimpsom(1/(x),1,3,3,1.0e-5) f =1.098620042680481step =8t =1.154407399999997图4在MATLAB中用变步长Simpson法计算j31 dx积分值1x2.2.3 Romberg 积分法Romberg 积分法是用 Richardson 外推算法来加快复合梯形求积公式的收敛速度，它的算 法如下，其中T是一系列逼近原定积分的Romberg积分值。m（1）计算b 一 aT（。）= f （a） + f （b）m2(2 j - 1)( b - a)2 k（2）对k = 1,2,3，计算下列各步：1b 一 a 2k-1T(k) =Tk-1 + 工 f (a +1212k 一1j=1对 m = 1,2,k 和 i = k, k 1, k 2,，计算:4mT i T i1i 1 =mmm +14 m 13）精度控制。面的计算过程如下表所示：表 1 Romberg 积分计算表格T(0)T(1)T (0)T (2)T (1)T (0)T (3)T (2)T (1)T (4)T (3)T (2)随着计算的步骤的增加，T( i)越来越逼近积分J bf (x) dx。下面是用T( m )来逼近 mmaJbf (x)dx 的 MATLAB 代码。a在MATLAB中编程实现的Romberg积分法的函数为：Romberg。功能： Romberg 积分法求函数的数值积分。调用格式： I,step=Romberg(f,a,b,eps) 。其中，f 为函数名；a 为积分下限；b 为积分上限；eps 为积分精度；I 为积分值；step 为积分划分的子区间次数。2.2.3.1 程序function I,step,t =Romberg(f,a,b,eps)%Romberg积分法求函数f在区间(a,b)上的定积分 %函数名： f%积分下限： a% 积分上限： b% 积分精度： eps%积分值： I%积分划分的子区间次数： stept0=cputime;if (nargin=3)eps=1.0e-4;end;M=1;tol=10;k=0;T=zeros(1,1);h=b-a;T(1,1)=(h/2)*(subs(sym(f),findsym(sym(f),a)+subs(sym(f),findsym(sym(f),b);% 初始值 while tolepsk=k+1;h=h/2;Q=0;for i=1:Mx=a+h*(2*i-1);Q=Q+subs(sym(f),findsym(sym(f),x);endT(k+1,1)=T(k,1)/2+h*Q;M=2*M;for j=1:kT(k+l,j+l)=T(k+l,j)+(T(k+l,j)-T(k,j)/(4Aj-l);endtol=abs(T(k+l,j+l)-T(k,j);endI=T(k+l,k+l);step=k;t=cputime-t0;%计算程序运行时间2232计算exdx积分值，精确到10-5：l在 MATLAB 命令窗口中输入下列命令： f,step,t=Romberg(exp(x),l,3,l.0e-5)f =l7.36725509503962lstep =4t =0.24960l5999999772.2.3.3计算fdx积分值，精确到10-5：lx在 MATLAB 命令窗口中输入下列命令： f,step,t=Romberg(l/(x),l,3,l.0e-5)f =l.0986l2289805927step =5t =0.436802800000009图5在MATLAB中用Romberg积分法计算j3 e”dx积分值 1图6在MATLAB中用Romberg积分法计算j 3丄dx积分值1x2.3 结果汇总由以上程序运行的结果列表如下：表 2 三种方法计算定积分对照表计算方法被求定积 分结果精度所用步 数用时（s）变步长梯形法exdxi17.367b031817855010-5106117.8431554000000变步长Simpson 法17.36726977894263991.528809800000005Romberg 积分 法17.36725509503962140.249601599999977变步长梯形法丄dx1 x1.098|9742247472310-54017.14450990000003变步长Simpson 法1.098|2004268048181.154407399999997Romberg 积分 法1.098|1228980592750.436802800000009注：这里面的步数指得是程序循环次数。第三章 结果分析3.1 结果简析由表2的汇总结果可以明显看出，在相同精度下（10-5）,以求积分J3exdx为例，变步1 长梯形法要计算106步，耗时117.8432秒，而且试过几次，用时差异并不十分明显。而变 步长Simpson法则只要9步，用时1.5288秒便得出了结果。Romberg积分法更是仅用4步， 用时0.2496秒便得出了结果。以上的时间是通过程序中调用cputime来计算时间的，虽然不 是十分精确（因为也包含了诸如选择、判断、赋值等时间，也即程序执行时间开始调用 cputime，来得到当时的时间，程序执行完毕时再调用cputime来求与程序开始时的时间差）， 但基本上只有ns数级的差异，并不影响结论的得出。确实，Romberg积分法（在这三种积 分法中）要更加具有速度上的优势。再说计算结果的精度问题，虽然精度都为 10-5，但结果上也存在着细微差异。与真实值 17.36725509472862 相比，变步长梯形的误差还要大些。表 3 三种方法计算定积分结果与真值对照表计算方法被求定 积分结果真实值精度所用 步数用时（s）变步长梯形法17.367|70318178550106117.843155400000017.3672155094728622505568242183229变步长Simpson 法J3 exdx117.36726977894263910-591.528809800000005Romberg 积分17.367255095040.2496015999法3962199977变步长梯形法1.0987974224747234017.144509900000031.0981228866变步长f3 丄 dx1 x1.0986200426810-581.1544073999Simpson 法04818109691395245236922599997Romberg 积分1.0986122898050.4368028000法592700009从表3三种方法计算定积分结果与真值对照中不难发现，变步长梯形法的计算精度是刚 好达到计算要求，而且其满足的精度是当前计算值与前一个计算值（程序中为I2-I1）的差 值小于10-5这个精度要求，从其与真实值的对比可以看出其精确度只达到了小数点后 3位， 而同样的前后计算数值的差值精度，变步长Simpson法则达到了小数点后4 位,而Romberg 积分法更是达到了小数点后8位之多！ Romberg积分法的优势非常明显，不仅是这三种数值 分析求定积分方法中在相同精度要求下与真实值相比精确位数最高的，而且也是最快的。因 此，就这三种方法中，如果精度要求较高的情况下， Romberg 积分法是非常理想的选择。3.2 原因探究那么究竟是什么原因造成这三种求定积分的方法有如此大的差异呢？首先，从原理与编程思想上着手来分析这个问题。将复合梯形求积公式中的每个子段x ,x 二等分，分点(新增节点)记为x ，得k k +1k + 0.5到2 n个子段x , x , x , x , x, x 。由复合梯形公式推出(步长为丄)0 0.50.51n - 0.5n2hA-1A-1T =一(f (a) + 2Z f (x ) + 2Z f (x ) + f (b)2 n 4kk + 0.5k 二 1k 二 0比较T与T得n2 n12n2nh n-1+ Z f (x )2k + 0.5k 二 0b-a(h =)n这个递推公式表明，求出T后，如果还要算T ，用上式计算只需补算新增分点x 上的 n2nk+0.5函数值 f (x) ，比直接用复合梯形求积公式计算计算量节约了一半。k + 0.5后误差公式，记h =，由n(0.5 h)2,I - T (f(b) - f(a)，2n12所以I T 1x ，1 -T 4n1I T x (T T )2 n 32n n根据上式，人们用二等分前后两次计算结果的差T T来估计T的误差。基本的计 2nn2 n算过程是先在a,b上用梯形公式算出T，然后把a,b二等分用递推公式算出T。计算二 12者误差，如果大于给定误差，继续计算T，检验T与T间的差来计算误差并决定是否继续3 3 2改变步长。变步长梯形法让计算机自动选择n的大小，使计算结果在允许的误差限内，这种算法属于自适应算法范畴，自适应算法的理论基础是后验误差估计。复合Simpson公式则是将子段x ,x 的中点记为x ，先在每个子段k k +1k + 0 . 5x , x (k = 0,1,n 1)上用Simpson公式计算积分k k +1xx xI k+1 f (x)dx - E IT(f (x ) + 4f (x _) + f (x )，x 6 k k+2 k+1然后把它们加起来并利用积分区间的可加性得I b f ( x ) dxahn 1n 1(f (a) + 4 Z f (x+ 2 Z f (x ) + f (b)=6k + 1k2k二 0k=l称 S 为复合 Simpson 积分值，也称上式为复合 Simpson 公式。n 下面比较收敛性与余项(截断误差)估计。先讨论复合梯形公式的收敛性。事实上，根 据积分的定义，复合梯形公式是用部分和作为定积分的一个近似值，所以当积分区间被分得很细(h很小)时，复合梯形公式求得的T有一定的精度，且当h T 0时，有T T I,即nn 用复合梯形公式求得的积分值收敛于准备积分值。复合梯形公式的余项(截断误差)I- T b f (x) dx = (f( b) f( a)； n12 a121h复合Simpson公式的余项(截断误差)I - S ()4(f(b) f(3)(a)。比较这n 180 2两个余项可以看到，当f (x)足够光滑时复合Simpson公式的精度比复合梯形公式高。这也就意味着，每次求出的积分值，复合Simpson公式都以比复合梯形公式更快的速度趋近真实 值，所以复合Simpson公式计算速度明显优于复合梯形公式。这一点由梯形公式与Simpson 公式(抛物线公式)的对比中也可见一斑！再来说 Romberg 积分法， Romberg 积分公式1C = R63 nn1C + 一2 n 63(C C )=2 nn64 C632 n其中C是复合4阶N-C积分值，n1161S +-(S S )=S S = C2 n 152 nn152 n15 nn其精度比 S 还要高。 还可以证明复合 4阶 N-C 公式的后验误差估计2nI C(C C )，用这个误差估计去修正C ，便得到了上面的Romberg积分公式。2 n 632 nn2 n可以证明R有且仅有7次代数精度，因此它属于机械求积公式的范畴但不是N-C型的。1所谓 Romberg 积分法是，在变步长过程中，按下图所示加工流程，用复合 Simpson 积分值、复合4阶N-C积分值及Romberg积分值把粗糙的梯形积分值T逐步加工成高精度的 nRomberg 积分值。1 TS21TSC421TSCR8421TSCR16 842逐行计算，算完前五行后得到R , R。若|R - R |W 8 （允许误差），就把R作为所求近似12 12 2值。否则再算第六行得R。直到前后两次Romberg值之差|r - R |w 8才停止计算，42 nn并把 R 作为所求近似值。2n由此可以， Romberg 积分法每一步的精度都要远高于变步长梯形法与变步长 Simpson 法，这就是为什么它仅用几步，也即仅计算几行便可以达到要求精度的原因。它的收敛速度 更加迅速。而且精度的提高不是一两位，而是好几位的大幅度提升。这也就是为什么在满足 的相同精度等级要求的同时， Romberg 积分法的精确度要明显优于变步长梯形法与变步长 Simpson 法。第四章 结论综合上述分析可以得出结论，当计算误差要求给定时，变步长Simpson法要明显快于变步长梯形法，而 Romberg 积分法则是精度要明显优于前两种算法，但速度上比变步长梯形法明显加快，但较之变步长Simpson法就没有比变步长梯形法那么明显，但还是有较大的提升。因此综合各种因素考虑，还是 Romberg 积分法有明显的速度及精度方面的优势。参考文献1.杨一都.数值计算方法.北京：高等教育出版社，2008.2.王正林，龚纯，何倩.精通 MATLAB 科学计算.北京：电子工业出版社，2007.