1、多种数值积分求积公式的分析比较吴春晖（中国海洋大学海洋环境学院 山东 青岛 266100）摘要：对于运用牛顿-莱布尼茨积分公式不能较好解决的定义在区间a,b上的可积函数，原函数并不能简单地用初等函数来表达，故需要构造定积分的近似计算公式。在本文中，主要构建了抛物线求积公式及其复化抛物线公式。在对抛物线类的求积公式进行应用检验后，再运用 Gauss 求积公式，构建 Gauss-Laguerre 求积公式，对相同的问题进行运用，并比较截断误差。之后再对求积过程进行优化，在限定误差范围的情况下，利用龙贝格算法，对求积加速收敛。关键词： 抛物线求积 复化求积 Gauss-Laguerre 加速收敛引言
2、：对于一些较为复杂的函数，在一定的误差要求下，需要通过构造的方式求给定函数的定积分。基本的替代法主要有梯形面积及抛物线近似代替曲边梯形。并通过划分更小的区间，减少截断误差从而提出了复化梯形及抛物线公式。为了提高运算效率，有加速收敛的Richardson 外推法和 Romberg 求积公式。之后，针对节点数固定情况下，提出了 Gauss 公式，使其获得最大的精度。本文主要研究的是抛物线求积法与 Gauss-Laguerre 公式。目录第一章 抛物线求积公式及应用 31.1 抛物线求积公式的算法 31.2 抛物线求积公式的 matlab 程序 31.3 复化抛物线求积公式的应用 .4第二章 Gau
3、ss-Laguerre 求积公式及应用 52.1 Gauss-Laguerre 的算法 .52.2Gauss-Laguerre 公式的 matlab 程序 .52.3 Gauss-Laguerre 求积公式的应用 .6第三章 龙贝格算法与算法优化 .73.1 龙贝格算法及程序 .73.2 利用龙贝格算法优化求积 .73.3 龙贝格算法的应用 8第四章 数值积分的分析总结 .9第一章 抛物线求积公式及应用1.1 抛物线求积公式的算法抛物线求积公式，是将区间二等分，以中点及两端点作为抛物线的三个点，并求出抛物线，在区间上对抛物线函数求积分。而复化抛物线公式，即在给定区间上，分成多个小份，并分别用抛
4、物线公式求积。所求得的数值积分余项如下：1.2 抛物线求积公式的 matlab 程序具体的程序代码如下：function s=simpr1(f,a,b,M) h=(b-a)/(2*M);s1=0;s2=0;for k=1:M x=a+h*(2*k-1);s1=s1+feval(f,x);end for k=1:(M-1)x=a+h*2*k; ,)(2180)( )(2180),(4 145bafhab xfhSdxfTfR knkkn 2,)(4)(2)(6)()(4)()()( 11102101 hxxfxfbfanbffxfhdfdf knknkkknkbak nkSfffdxf 1)()
5、()()(s2=s2+feval(f,x);end s=h*(feval(f,a)+feval(f,b)+4*s1+2*s2)/3; 说明：f 代表的是原函数，该程序为抛物线复化积分公式。当 M 值为 1 时，则为简化的抛物线求积公式。最后 s 即为用 matlab 语言所表述的抛物线求积公式。1.3 复化抛物线求积公式的应用在应用部分，所计算的积分为 1/(1+x2)。将函数与节点，积分区间输入函数后，得出的结果如下图表所示。图 1.3.1 复化抛物线求积（所取节点为 5-50）Fig.1.3.1Compound parabolic quadrature (from the node for
6、 5-50)图 1.3.2 积分误差（节点数为 0-50）Fig.1.3.2quadrature error (from the node for 0-50)节点数0 5 10 15 20 25 30 35 40 45 50误差 0.0020659.91E-091.55E-101.36E-112.42E-126.35E-132.13E-138.44E-143.77E-141.87E-141.01E-14表 1.3.1 每隔 5 节点的误差Table 1.3.1 every 5 node error分析总结：复化抛物线积分对给定的函数有较好的适用性，在节点数为 3 时，误差就处于 E-3 数量级
7、。这主要与函数的特性及选择区间有关，所求得积分为 1/4 圆，故抛物线较能吻合特征。从图像我们可以观察到：随着节点数的增加，积分的准确性得到增强，且增加速度较快。在 50 个节点的情况下，误差达到了 E-14 的数量级。从误差观察得出，误差减少到下一数量级的节点间隔数增大。第二章 Gauss-Laguerre 求积公式及应用2.1 Gauss-Laguerre 的算法Gauss-Laguerre 求积公式是 Gauss 求积公式的一种建立在无穷区间上的特殊求积公式。Laguerre 多项式：在 关于权函数为 的正交多项式0,)()xe（1 ） 。()nxxndeL故在求积分时，我们主要使用的求
8、积公式为： （2）01()()nxkefdAfx其中, （k=1,n）是 的 n 个零点，求积系数 k=1,2,n （3）kx()Lx2!()knkL2.2Gauss-Laguerre 公式的 matlab 程序程序代码如下：syms xf;syms l;n=20;l(n)=exp(x)*diff(x(n)*exp(-x),x,n);xs=solve(char(l(n),=0,x);ll=diff(l(n),1);for i=1:n a(i)=(factorial(n)2/(double(xs(i)*(subs(ll,x,double(xs(i)2);ends=double(a)*feval(
9、f1,double(xs);程序分析：n代表Laguerre的项数，f 为函数。该程序先求出相应的Laguerre多项式后，通过solve 解出节点，再用节点算出系数，代入求积公式。2.3 Gauss-Laguerre 求积公式的应用所选择的积分所在的区间为0 inf，f(x)= x，取不同的项数，进行求积，将所得的结果整理如下。1 2 3 5 10 150.92388 0.90644 0.89928 0.893296 0.889062 0.88784120 25 28 29 30 0.887299 0.887005 0.886962 0.886923 0.886887 表 2.3.1 积分数
10、据表Table.2.3.1 the data table图 2.2.2 Gauss-Laguerre 公式积分曲线Figure.2.2.2 Gauss-Laguerrequadrature分析总结：Gauss-Laguerre 公式的适用对象为在特定区间，即0,inf区间上的权函数乘以f(x)形式的特定积分。从图像可以观察得到，Gauss-Laguerre 公式随着节点数的增大，呈双曲线的形状。其积分精度逐渐提高，但收敛较慢。而用其它形式的积分来处理无限区间的积分，积分的计算量更大，精度也比较差。所以 Gauss-Laguerre 公式对特定积分的效果较好，但要达到一定精度，计算节点数需要到一
11、个较大值。第三章 龙贝格算法与算法优化3.1 龙贝格算法及程序龙贝格算法从简单的梯形序列开始逐步进行线性加速，具有占用内存少，精度高的优点。而梯形或抛物线的积分算法的区间并不是线性增加的，每次计算的效率会越来越低。故龙贝格算法适用于给定精度下的较为高效的求积。龙贝格的 matlab 程序如下：fuction t=Romberg(fun,a,b,e)if narginei=i+1;h=h/2;T(i+1,1)=T(i,1)/2+sum(feval(fun,a+h/2:h:b-h/2+0.001*h)*h/2;for j=1:i;T(i+1,j+1)=4j*T(i+1,j)/(4j-1)-T(i,
12、j)/(4j-1);endendt=T(i+1,j+1);3.2 利用龙贝格算法优化求积在之前的内容中，利用抛物线复化公式计算了 1./(1+x.2)，在区间1,2 上的积分为 pi/4,在不同的节点取值下，抛物线复化求积的误差精度不同。如果要得到 e-10 的精度范围，需要取的节点数为 8。而选用龙贝格算法能有效地减少计算步骤，提高效率。选取精度为 e-10e-16，分别计算对应的复化求积公式和龙贝格法算法所需要运行的步骤。选取函数为 1./(1+x.2)，两种算法的结果如下：精度 10 11 12 13 14 15 16复化求积 8 11 16 24 35 51 72龙贝格算法 7 8 8
13、 8 9 9 11表 3.2.1 计算次数与求积公式的关系Table 3.2.1 The number calculation, and the quadrature formula图 3.2.2 计算次数与求积公式的关系Figure.3.2.2 The number calculation, and the quadrature formula从图像上，可以发现龙贝格的算法效率比抛物线的复化求积公式高，且随着精度要求的提高，抛物线的复化求积所取得等分点增加得越来越快。而龙贝格算法的次数平滑增长，在较小的次数，就能达到极高的精度。3.3 龙贝格算法的应用龙贝格算法对限定区间的求积较为适用。在给
14、定数据精度的情况下，对 f(x)=sqrt(x)在1,9的区间上积分。所要求的精度为 10-5.将代数值输入函数中，计算所得到的值为17.3333，符合要求。随着精度的提升，应用龙贝格算法。对给定函数在特定区间上的积分的值绘成的曲线如下。图 3.3.1 Romberg 算法积分曲线Figure.3.3.1 Romberg integral curve algorithm从图像中我们可以发现，由于设置了精度，曲线的攀升效果较为明显，且较大精度值精度过高，基本没有区别。故龙贝格算法具有较好的准确性与较高的效率。第四章 数值积分的分析总结这篇文章主要对几种基本求积公式进行了 matlab 的程序实现
15、与相应的应用，并简单地分析对比了各个求积公式的优缺点，归纳如下：复化梯形或抛物线求积公式的应用较为广泛，实现手段较为简单。在对精度要求不大，积分区间为一小段实数的情况下，能较好地实现求积。求积的结果能够粗略地反映积分的大小。且随着节点数的增加，数值积分的值逼近原值。在该文中，主要用的是抛物线的复化积分，对曲线的求积效果可能较好，而对折线段的求积可能没有梯形求积公式精准。Gauss 求积的计算过程比较繁琐，在本文中的 Gauss-Laguerre 公式求积公式需要先计算节点与对应系数，在取这些点的值，代入函数来积分。运算较多，可能需要在求积前，需要预先导入节点的值与相应的系数来优化过程。Gaus
16、s 的求积公式可以运用在范围较大的或者无限的区间。而其余的求积公式并不能很好实现在一个较大区间或特定区间的取值。Gauss 公式的数据精度较高。龙贝格算法的优点上文有多次提及，即通过较为高效的方式求出精度较高的积分。在通过与复化抛物线积分的比较时，我们发现龙贝格在求得指定精度值，所需要进行的计算次数，远小于抛物线复化求积公式所要求的划分节点数，且随着精度要求的提高，增长较为平缓。在具体的数值积分情形下，可以利用多种求积方式来解决问题。如基于抛物线积分的自适应抛物线法。其主要原理是对子区间进行了分别积分。定下误差份额，若部分满足，则保留。对不满足的部分，切分为较小区间再积分，值得所有子区间均符合
17、积分条件。自适应的数值积分还可以选用其他的积分，如龙贝格算法。但这样某种程度上增加了计算难度。在对精度要求极为严苛时，误差判断是必须的，而如何确定误差，或者说是余项，各个积分公式都有对应的截断误差。对截断误差随 n 值的变化，可以发现。Gauss 公式及龙贝格法的收敛快于复化抛物线。对于数值积分，选取好的公式与算法的优化较为重要。参考文献1 徐萃薇，孙绳武.计算方法引论.北京：高等教育出版社，20152 郑成德。数值计算方法.北京：清华大学出版社，2010Abstract:For the use of Newton integral formula cant better solve defi
18、ned on the interval (a, b) integrable function, the function is not simply to express with elementary function, so the need to construct an approximate calculation formula of definite integral In this paper, the main building of parabolic quadrature formula and its compound parabolic equation Applic
19、ation in quadrature formula of parabolic type, after testing and using the Gauss quadrature formula, build the Gauss - Laguerre quadrature formula, to use the same problem, and compare the truncation error and then optimize the quadrature process, in the case of a limited range of error, the quadrature accelerating convergence Keywords: parabolic quadrature After the quadrature Gauss Laguerre accelerating convergence
