多种数值积分求积公式的分析比较
吴春晖
（中国海洋大学海洋环境学院山东青岛 266100）
摘要：
对于运用牛顿-莱布尼茨积分公式不能较好解决的定义在区间[a,b]上的可积函数，原函数并不能简单地用初等函数来表达，故需要构造定积分的近似计算公式。在本文中，主要构建了抛物线求积公式及其复化抛物线公式。在对抛物线类的求积公式进行应用检验后，再运用Gauss求积公式，构建Gauss-Laguerre求积公式，对相同的问题进行运用，并比较截断误差。之后再对求积过程进行优化，在限定误差范围的情况下，利用龙贝格算法，对求积加速收敛。
关键词：抛物线求积复化求积Gauss-Laguerre加速收敛
引言：
对于一些较为复杂的函数，在一定的误差要求下，需要通过构造的方式求给定函数的定积分。基本的替代法主要有梯形面积及抛物线近似代替曲边梯形。并通过划分更小的区间，减少截断误差从而提出了复化梯形及抛物线公式。为了提高运算效率，有加速收敛的Richardson外推法和Romberg求积公式。之后，针对节点数固定情况下，提出了Gauss公式，使其获得最大的精度。本文主要研究的是抛物线求积法与Gauss-Laguerre公式。
目录
第一章抛物线求积公式及应用 (3)
1.1抛物线求积公式的算法 (3)
1.2抛物线求积公式的matlab程序 (3)
1.3复化抛物线求积公式的应用 (4)
第二章Gauss-Laguerre求积公式及应用 (5)
2.1 Gauss-Laguerre的算法 (5)
2.2Gauss-Laguerre公式的matlab程序 (5)
2.3 Gauss-Laguerre求积公式的应用 (6)
第三章龙贝格算法与算法优化 (7)
3.1龙贝格算法及程序 (7)
3.2利用龙贝格算法优化求积 (7)
3.3 龙贝格算法的应用 (8)
第四章数值积分的分析总结 (9)
第一章 抛物线求积公式及应用
1.1 抛物线求积公式的算法
抛物线求积公式，是将区间二等分，以中点及两端点作为抛物线的三个点，并求出抛
物线，在区间上对抛物线函数求积分。而复化抛物线公式，即在给定区间上，分成多个小份，并分别用抛物线公式求积。
所求得的数值积分余项如下：
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;
]
,[)(2
180)()
(2
180)(),()4(441
10)4(4
5
b a f h a b x x f h S dx x f T f R k k n k k b a n n ∈?--=≤≤?-=-=+-=∑
?
ηηηη2,))(4)(2)()((6))
()(4)((6)()(2
11110
2
112
11
0101
h
x x
x f x f b f a f n a b x f x f x f h dx x f dx x f k k n k n k k k k k k n k b a n k x x k k
+
=+++-=++≈=+
-=-=+++-=-=∑∑∑
?
∑?
+n
n k n k k k b a S x f x f b f a f n
a b dx x f =+++-≈∑∑?
-=-=+1
1
1
2
1))(4)(2)()((6)(
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+x^2)。将函数与节点，积分区间输入函数后，得出的结果如下图表所示。
图1.3.1 复化抛物线求积（所取节点为5-50）
Fig.1.3.1Compound parabolic quadrature (from the node for 5-50)
图1.3.2 积分误差（节点数为0-50）
Fig.1.3.2quadrature error (from the node for 0-50)
Table 1.3.1 every 5 node error
分析总结：复化抛物线积分对给定的函数有较好的适用性，在节点数为3时，误差就处于E-3数量级。这主要与函数的特性及选择区间有关，所求得积分为1/4圆，故抛物线较能吻合特征。从图像我们可以观察到：随着节点数的增加，积分的准确性得到增强，且增加速度较快。在50个节点的情况下，误差达到了E-14的数量级。从误差观察得出，误差减少到下一数量级的节点间隔数增大。
第二章Gauss-Laguerre 求积公式及应用
2.1 Gauss-Laguerre 的算法
Gauss-Laguerre 求积公式是Gauss 求积公式的一种建立在无穷区间上的特殊求积公式。Laguerre 多项式：在[0,)-∞关于权函数为()x
x e ρ-=的正交多项式
()
()n n x x
n n
d x
e L x e dx -=（1）。
故在求积分时，我们主要使用的求积公式为：
1
()()n
x
k k k e f x dx A f x ∞
-=≈∑?
（2）
其中,k x （k=1,…,n ）是()n L x 的n 个零点，求积系数2
'2
(!)[()]k n k k
n A L x x =k=1,2,…,n （3）
2.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) );
end
s=double(a)*feval(f1,double(xs));
程序分析：n代表Laguerre的项数，f为函数。该程序先求出相应的Laguerre多项式后，通过solve解出节点，再用节点算出系数，代入求积公式。
2.3 Gauss-Laguerre求积公式的应用
所选择的积分所在的区间为[0 inf]，f(x)= √x，取不同的项数，进行求积，将所得的结果整理如下。
表2.3.1 积分数据表
Table.2.3.1 the data table
图2.2.2Gauss-Laguerre公式积分曲线
Figure.2.2.2Gauss-Laguerrequadrature
分析总结：Gauss-Laguerre公式的适用对象为在特定区间，即[0,inf]区间上的权函数乘以f(x)形式的特定积分。从图像可以观察得到，Gauss-Laguerre公式随着节点数的增大，呈双曲线的形状。其积分精度逐渐提高，但收敛较慢。而用其它形式的积分来处理无限区间的积分，积分的计算量更大，精度也比较差。所以Gauss-Laguerre公式对特定积分的效果较好，但要达到一定精度，计算节点数需要到一个较大值。
第三章龙贝格算法与算法优化
3.1龙贝格算法及程序
龙贝格算法从简单的梯形序列开始逐步进行线性加速，具有占用内存少，精度高的优点。而梯形或抛物线的积分算法的区间并不是线性增加的，每次计算的效率会越来越低。故龙贝格算法适用于给定精度下的较为高效的求积。
龙贝格的matlab程序如下：
fuction t=Romberg(fun,a,b,e)
if nargin<4
e=1e-8;
end
i=1;j=1;
h=b-a;
T(i,1)=h/2*(feval(fun,a)+feva1(fun,b));
T(i+1,1)=T(i,1)/2+sum(feval(fun,a+h/2:h:b-h/2+0.001*h))*h/2;
T(i+1,j+1)=4.^j.*T(i+1,j)/(4.^j-1)-T(i,j)/(4^j-1);
while abs(T(i+1,i+1)-T(i,i))>e
i=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)=4^j*T(i+1,j)/(4^j-1)-T(i,j)/(4^j-1);
end
end
t=T(i+1,j+1);
3.2利用龙贝格算法优化求积
在之前的内容中，利用抛物线复化公式计算了1./(1+x.^2)，在区间[1,2]上的积分为pi/4,在不同的节点取值下，抛物线复化求积的误差精度不同。如果要得到e-10的精度范围，需要取的节点数为8。而选用龙贝格算法能有效地减少计算步骤，提高效率。
选取精度为e-10~e-16，分别计算对应的复化求积公式和龙贝格法算法所需要运行的步骤。
选取函数为1./(1+x.^2)，两种算法的结果如下：
表3.2.1计算次数与求积公式的关系
Table 3.2.1 The number calculation, and the quadrature formula
图3.2.2计算次数与求积公式的关系
Figure.3.2.2The number calculation, and the quadrature formula
从图像上，可以发现龙贝格的算法效率比抛物线的复化求积公式高，且随着精度要求的提高，抛物线的复化求积所取得等分点增加得越来越快。而龙贝格算法的次数平滑增长，在较小的次数，就能达到极高的精度。
3.3 龙贝格算法的应用
龙贝格算法对限定区间的求积较为适用。在给定数据精度的情况下，对f(x)=sqrt(x)在[1,9]的区间上积分。所要求的精度为10^-5.将代数值输入函数中，计算所得到的值为17.3333，
符合要求。
随着精度的提升，应用龙贝格算法。对给定函数在特定区间上的积分的值绘成的曲线如下。
图3.3.1Romberg算法积分曲线
Figure.3.3.1 Romberg integral curve algorithm
从图像中我们可以发现，由于设置了精度，曲线的攀升效果较为明显，且较大精度值精度过高，基本没有区别。故龙贝格算法具有较好的准确性与较高的效率。
第四章数值积分的分析总结
这篇文章主要对几种基本求积公式进行了matlab的程序实现与相应的应用，并简单地分析对比了各个求积公式的优缺点，归纳如下：
复化梯形或抛物线求积公式的应用较为广泛，实现手段较为简单。在对精度要求不大，积分区间为一小段实数的情况下，能较好地实现求积。求积的结果能够粗略地反映积分的大小。且随着节点数的增加，数值积分的值逼近原值。在该文中，主要用的是抛物线的复化积分，对曲线的求积效果可能较好，而对折线段的求积可能没有梯形求积公式精准。
Gauss求积的计算过程比较繁琐，在本文中的Gauss-Laguerre公式求积公式需要先计算节点与对应系数，在取这些点的值，代入函数来积分。运算较多，可能需要在求积前，需要预先导入节点的值与相应的系数来优化过程。Gauss的求积公式可以运用在范围较大的或者无限的区间。而其余的求积公式并不能很好实现在一个较大区间或特定区间的取值。Gauss 公式的数据精度较高。
龙贝格算法的优点上文有多次提及，即通过较为高效的方式求出精度较高的积分。在通过与复化抛物线积分的比较时，我们发现龙贝格在求得指定精度值，所需要进行的计算次数，远小于抛物线复化求积公式所要求的划分节点数，且随着精度要求的提高，增长较为平缓。
在具体的数值积分情形下，可以利用多种求积方式来解决问题。如基于抛物线积分的自
适应抛物线法。其主要原理是对子区间进行了分别积分。定下误差份额，若部分满足，则保留。对不满足的部分，切分为较小区间再积分，值得所有子区间均符合积分条件。
自适应的数值积分还可以选用其他的积分，如龙贝格算法。但这样某种程度上增加了计算难度。
在对精度要求极为严苛时，误差判断是必须的，而如何确定误差，或者说是余项，各个积分公式都有对应的截断误差。对截断误差随n值的变化，可以发现。Gauss公式及龙贝格法的收敛快于复化抛物线。
对于数值积分，选取好的公式与算法的优化较为重要。
参考文献
[1] 徐萃薇，孙绳武.计算方法引论.北京：高等教育出版社，2015
[2] 郑成德。数值计算方法.北京：清华大学出版社，2010
Abstract:
For the use of Newton integral formula can't better solve defined 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 Application 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
单元节点和积分点有什么区别 学过数值积分的应该知道，有限元中的积分点指高斯积分点，因为这些点的收敛性好，精度高。 1、节点 在单元内，采用形函数来表述单元内变量的分布规律。而节点值是在节点处的对应物理量。 以简单矩形单元的温度为例： 四个节点i,j,m,n的温度分别为Ti,Tj,Tm,Tn. 则以单元内自然坐标(x,y),(-1,-1),(-1,1),(1,-1),(1,1)分别为四个节点，单元内温度分布为： T={Si, Sj, Sm, Sn} {Ti, Tj, Tm, Tn} Si=1/4(1-x)(1-y) Sj=1/4(1+x)(1-y) Sm=1/4(1+x)(1+y) Sn=1/4(1-x)(1+y) (单元的形函数我们可以从手册中查到) 从而我们知道了温度在单元内的分布。 2、积分节点 我们需要对温度在单元内的面积上进行积分时，因为节点的温度显然与x,y无关，我们只需要考虑对形函数积分。 采用Gauss_Legendre多项式计算积分时，我们只需要计算根据特定积分点的值（在自然坐标系下是固定的，可以查手册，这些点也叫高斯点、积分点）并加以权重就可以。这就把复杂的积分问题变成了简单的代数问题。因为形函数只有单元有关，所以积分点也只与单元形状有关。 3.应力一般采用多个积分点的相互插值或外延来计算节点应力。这只是为了减少误差。因为在积分点应力比节点具有更高阶的误差。 从理论上说，形函数已知后，用Maple或者Mathematic等软件进行符号积分的话，是可以精确计算出刚度矩阵和质量矩阵，但是这样做的话，对于工程实际应用来说并不合适（ 原因：1，费时；2，Mindlin中厚板有剪力锁死问题，有时候需要采用缩聚积分），所以有些书上会把2节点梁单元的刚度阵直接写出来，但是再复杂点的单元，就使用数值积分（Newton-Cotes积分和高斯积分） 高斯积分的话，积分点不在节点上 牛顿-科斯的积分点就是节点，这样得到的质量矩阵是集中质量阵形式 个人理解： 1.节点作用构造形函数，节点的多少描述规则形状单元内的应力的近似分布情况，并获取节点上的位移值
高斯—勒让德积分公式 摘要： 高斯—勒让德积分公式可以用较少节点数得到高精度的计算结果,是现在现实生活中经常运用到的数值积分法。然而，当积分区间较大时，积分精度并不理想。 T he adva ntage of Gauss-Legendre integral formula is tend to get high-precision calculational result by using fewer Gauss-points, real life is now often applied numerical integration method. But the precision is not good when the length of integral interval is longer. 关键字： … 积分计算，积分公式，高斯—勒让德积分公式，MATLAB Keyword: Integral Calculation , Integral formula ,Gauss-Legendre integral formula, Matlab 引言： 众所周知，微积分的两大部分是微分与积分。微分实际上是求一函数的导数，而积分是已知一函数的导数，求这一函数。所以，微分与积分互为逆运算。 】 实际上，积分还可以分为两部分。第一种，是单纯的积分，也就是已知导数求原函数，称为不定积分。 相对而言，另一种就是定积分了，之所以称其为定积分，是因为它积分后得出的值是确定的，是一个数，而不是一个函数。 计算定积分的方法很多，而高斯—勒让德公式就是其中之一。 高斯积分法是精度最高的插值型数值积分，具有2n+1阶精度，并且高斯积分总是稳定。而高斯求积系数，可以由Lagrange多项式插值系数进行积分得到。 高斯—勒让德求积公式是构造高精度差值积分的最好方法之一。他是通过让节点和积分系数待定让函数f(x)以此取i=0,1,2....n次多项式使其尽可能多的能够精确成立来求出积分节点和积分系数。高斯积分的代数精度是2n-1，而且是最
高斯型积分公式-CAL-FENGHAI.-(YICAI)-Company One1
2 Guass-Legendre 积分程序 1. 目的意义： 可以提高数值积分的代数精度 2. 数学公式： ) ()()(1k n k k b a x f A dx x f x ∑?=≈ρ 3. 程序： #include<> #include<> #define N 10 #define f(x) (cos(x)) int main() { int n=0; int k=0; int i=0; double x[N]={}; double A[N]={}; double s=; n=2; switch(n)
{ case 1: { x[1]=0; A[1]=2; break; } case 2: { x[1]=; x[2]=; A[1]=1; A[2]=1; break; } case 3: { x[1]=; x[2]=0; x[3]=; A[1]=; A[2]=; 3
A[3]=; break; } case 4: { x[1]=; x[2]=; x[3]=; x[4]=; A[1]=; A[2]=; A[3]=; A[4]=; break; } default: { printf("error! 请添加数据！\n"); return 0; } } 4
for(i=1;i<=n;i++) { s=s+A[i]*f(x[i]); } printf("由高斯-勒让德积分公式计算得I=%lf\n",s); return 0; } 4.运行结果： 5.参考文献: [1] 谭浩强. C语言程序设计[M]. 北京：清华大学出版社，2005. [2] 秦新强. 数值逼近, 西安,2010. 5
摘要 求函数在给定区间上的定积分，在微积分学中已给出了许多计算方法，但是，在实际问题计算中，往往仅给出函数在一些离散点的值，它的解析表达式没有明显的给出，或者，虽然给出解析表达式，但却很难求得其原函数。这时我们可以通过数值方法求出函数积分的近似值。 当然再用近似值代替真实值时，误差精度是我们需要考虑因素，但是除了误差精度以外，还可以用代数精度来判断其精度的高低。已知n+1点的Newton-Cotes型积分公式，当n为奇数时，其代数精度为n；当n为偶数时，其代数精度达到n+1。若对随机选取的n+1个节点作插值型积分公式也仅有n次代数精度。 如何选取适当的节点，能使代数精度提高？Gauss型积分公式可是实现这一点，但是Gauss型求积公式，需要被积函数满足的条件是正交，这一条件比较苛刻。因此本实验将针对三种常用的Gauss型积分公式进行讨论并编程实现。 关键词：Newton-Cotes型积分公式正交多项式代数精度
1、实验目的 1)通过本次实验体会并学习Gauss型积分公式，在解决如何取节点能提 高代数精度这一问题中的思想方法。 2)通过对Gauss型积分公式的三种常见类型进行编程实现，提高自己的 编程能力。 3)用实验报告的形式展现，提高自己在写论文方面的能力。 2、算法流程 下面介绍三种常见的Gauss型积分公式 1)高斯-勒让德（Gauss-Legendre）积分公式 勒让德（Legendre）多项式 如下定义的多项式 称作勒让德多项式。由于是次多项式，所以是n次多项式，其最高次幂的系数与多项式 的系数相同。也就是说n次勒让德多项式具有正交性即勒让德多项式 是在上带的n次正交多项式，而且 这时Gauss型积分公式的节点就取为上述多项式的零点，相应的Gauss型积分公式为 此积分公式即成为高斯-勒让德积分公式。 其中Gauss-Legendre求积公式的系数
实习论文 题目高斯勒让德积分公式 专业信息与计算科学 班级计算092 学号3090811065 学生周吉瑞 指导教师秦新强 2011 年
高斯勒让德积分公式 专 业： 信息与计算科学 学 生： 周吉瑞 指导老师： 秦新强 摘要 关于数值积分公式0 ()()b n k k k a f x dx A f x =≈∑?，除了用误差来分析其精度以外，还可以 用代数精度来判断其代数精度的高低，已知n+1点Newton-Cotes 型积分公式，当n 为奇数时，其代数精度为n ，当n 为偶数时，其代数精度达到n+1。 n+1点的Newton-Cotes 型积分公式属于插值积分型积分公式，一般地，若对随机选取的n+1个节点作插值型积分公式也仅有n 次代数精度，但是，如果求积节点选取适当，就有可能提高数值积分的代数精度，高斯型积分公式就可以实现这一目标。 关 键 词：数值积分，代数精度，高斯型积分公式
一、目的意义 构造Gaoss 型求积公式除需要求出正交多项式外，还需要求出正交多项式的零点和求积系数，当3n ≥时，这些工作均很困难，因此给出高斯-勒让德积分公式的零点和系数。 二、公式 高斯-勒让德积分公式：1 1 1 ()()n k k k f x A f x -=≈∑?； 三、算法流程 Step1:输入所用的点数n ； Step2:对i=1，2，···，n 循环执行步3； Step3:I= I+ ()i i A f x ; Step4:输出I ；结束。
四、算法程序 #include
Gauss型积分公式
摘要 求函数在给定区间上的定积分，在微积分学中已给出了许多计算方法，但是，在实际问题计算中，往往仅给出函数在一些离散点的值，它的解析表达式没有明显的给出，或者，虽然给出解析表达式，但却很难求得其原函数。这时我们可以通过数值方法求出函数积分的近似值。 当然再用近似值代替真实值时，误差精度是我们需要考虑因素，但是除了误差精度以外，还可以用代数精度来判断其精度的高低。已知n+1点的Newton-Cotes型积分公式，当n为奇数时，其代数精度为n；当n 为偶数时，其代数精度达到n+1。若对随机选取的n+1个节点作插值型积分公式也仅有n次代数精度。 如何选取适当的节点，能使代数精度提高？Gauss型积分公式可是实现这一点，但是Gauss型求积公式，需要被积函数满足的条件是正交，这一条件比较苛刻。因此本实验将针对三种常用的Gauss型积分公式进行讨论并编程实现。 关键词：Newton-Cotes型积分公式正交多项式代数精度
1、实验目的 1)通过本次实验体会并学习Gauss型积分公式，在解决如何取节点能提 高代数精度这一问题中的思想方法。 2)通过对Gauss型积分公式的三种常见类型进行编程实现，提高自己的 编程能力。 3)用实验报告的形式展现，提高自己在写论文方面的能力。 2、算法流程 下面介绍三种常见的Gauss型积分公式 1)高斯-勒让德（Gauss-Legendre）积分公式 勒让德（Legendre）多项式 如下定义的多项式 称作勒让德多项式。由于是次多项式，所以是n次多项式，其最高次幂的系数与多项式 的系数相同。也就是说n次勒让德多项式具有正交性即勒让德多项式是在上带的n次正交多项式，而且 这时Gauss型积分公式的节点就取为上述多项式的零点，相应的Gauss型积分公式为 此积分公式即成为高斯-勒让德积分公式。 其中Gauss-Legendre求积公式的系数 1
实习论文 题目高斯拉盖尔积分公式 专业信息与计算科学 班级计算092 学号3090811065 学生周吉瑞 指导教师秦新强 2011 年
高斯拉盖尔积分公式 专 业： 信息与计算科学 学 生： 周吉瑞 指导老师： 秦新强 摘要 关于数值积分公式0()()b n k k k a f x dx A f x =≈∑?，除了用误差来分析其精度以外，还可以 用代数精度来判断其代数精度的高低，已知n+1点Newton-Cotes 型积分公式，当n 为奇数时，其代数精度为n ，当n 为偶数时，其代数精度达到n+1。 n+1点的Newton-Cotes 型积分公式属于插值积分型积分公式，一般地，若对随机选取的n+1个节点作插值型积分公式也仅有n 次代数精度，但是，如果求积节点选取适当，就有可能提高数值积分的代数精度，高斯型积分公式就可以实现这一目标。 关 键 词：数值积分，代数精度，高斯型积分公式
一、目的意义 构造Gaoss 型求积公式除需要求出正交多项式外，还需要求出正交多项式的零点和求积系数，当3n ≥时，这些工作均很困难，因此给出高斯-拉盖尔积分公式的零点和系数。 二、公式 高斯-拉盖尔积分公式：10()()n x k k k e f x A f x ∞-=≈∑?； 三、算法流程 Step1:输入所用的点数n ； Step2:对i=1，2，···，n 循环执行步3； Step3:I= I+ ()i i A f x ; Step4:输出I ；结束。 四、算法程序 #include
数学软件实验任务书
实验1Romberg 积分法 1实验原理 Romberg 方法是实用性很强的一种数值积分方法，其收敛速度是很快的，这里给出Romberg 积分的计算方法。 （1）计算1 (0,0)()[()()]2 R b a f a f b =-+ （2）计算2 211111 (,0)(1,0)(())222 i i i k h R i i f a k h ---==-++-∑ （3）计算114(,1)(1,1)(,)41 j j R m j R m j R m j ------=- 2 实验数据 用 Romberg 积分方法计算： 15 2 4x S dx x =+? 3实验程序 程序1 function s=rombg(a,b,TOL) n=1; h=b-a; delt=1; x=a; k=0; R=zeros(4,4); R(1,1)=h*(rombg_f(a)+rombg_f(b))/2; while delt>TOL k=k+1; h=h/2; s=0; for j=1:n x=a+h*(2*j-1); s=s+rombg_f(x); end R(k+1,1)= R(k,1)/2+h*s; n=2*n; for i=1:k R(k+1,i+1)=((4^i)*R(k+1,i)-R(k,i))/(4^i-1); end delt=abs(R(k+1,k)-R(k+1,k+1)); end
s=R(k+1,k+1); 程序2 function f=rombg_f(x) f=x/(4+x^2); 程序3 s=rombg(0,1.5,1.e-6) %作出图形 x=0:0.02:1.5; y=x./(4+x.^2); area(x,y) grid 4 实验结果 s = 0.2231 实验2高斯-勒让德积分法1实验原理 Gauss-Legendre求积公式为
§4-4 高斯积分法及其应用 ● 由§4-3知，在计算空间等参数单元的载荷列阵及刚度矩阵时，需用到如下形式的定积分： ηξηξd d f ? ? --1 11 1 ),(； ζηξζηξd d d f ? ?? ---1 1111 1 ),,( 其中被积分函数f(ξ，η，ζ)一般是很复杂的，即使能够得出它的显式，其积分也是很繁的。因此，一般用数值积分来代替函数的定积分。 ● 数值积分：在积分区域内按一定规则选出一些点，称为积分点，算出被积函数f(ξ， η，ζ)在这些积分点处的值，然后再乘以相应的加权系数并求和，作为近似的积分值。 ● 数值积分的方法有多种，其中高斯积分法可以用相同的积分点数达到较高的精度，或者说用较少的积分数达到同样的精度。 一、高斯积分法 1．一维积分的高斯公式 ● 一维积分的高斯公式 ∑? =-= n i i i f H d f 1 1 1 )()(ξξξ (4-47) 其中f(ξi )是被积函数在积分点ξi 处的数值，H i 为加数系数，n 为积分点数目。 ● 可以证明， ? 对于n 个积分点，只要选取适当的加数系数及积分点位置，能够使(4-47)式在被积分函数为不超过(2n-1)次多项式时精确成立。 ? 由于多数函数可表示成多项式形式，这种积分适应于大多数函数。 ● 例如， ? n=1时 )()(111 1 ξξξf H d f I == ? - (a) 不论f(ξ)的次数是0还是1，只需取H １=2，ξ１ ＝０，上式均是精确成立的。因为 ξξ10)(C C f += (b)
101 ()22(0)I f d C f ξξ-= ==?? (c) ? 当n=2时，能保证(4-47)式精确成立所允许的多项式的最高次数是3，此时，f(ξ)的通式为 3 32 210)(ξξ ξξC C C C f +++= (d) 其精确积分为 201 1 3 22)(C C d f I + == ? -ξξ (e) 数值积分为 ) ()() ()()(3 232 221023 132 1210122112 1 ξξξξξξξξξC C C C H C C C C H f H f H f H I i i i +++++++=+== ∑= (f) 为了在C 0～C 3取任意值(包括取零值在内)时公式(f)是精确的，显然应有 221=+H H ， 02211=+ξξH H 3 22 222 11= +ξξH H ， 0322311=+ξξH H 所以，应取 2,269,350,577.03 121-=- =-=ξξ 0,000,000,000.121==H H ? 同样，对于不超过五次的多项式，只要取 n=3 130.577,350,269,2ξξ=-=- =- 0,000,000,000.02=ξ 6,555,555,555.09 521== =H H 9,888,888,888.09 83== H
Guass-Legendre 积分程序 1. 目的意义： 可以提高数值积分的代数精度 2. 数学公式： ) ()()(1k n k k b a x f A dx x f x ∑?=≈ρ 3. 程序： #include<> #include<> #define N 10 #define f(x) (cos(x)) int main() { int n=0; int k=0; int i=0; double x[N]={}; double A[N]={}; double s=; n=2; switch(n)
{ case 1: { x[1]=0; A[1]=2; break; } case 2: { x[1]=; x[2]=; A[1]=1; A[2]=1; break; } case 3: { x[1]=; x[2]=0; x[3]=; A[1]=; A[2]=;
A[3]=; break; } case 4: { x[1]=; x[2]=; x[3]=; x[4]=; A[1]=; A[2]=; A[3]=; A[4]=; break; } default: { printf("error! 请添加数据！\n"); return 0; } }
for(i=1;i<=n;i++) { s=s+A[i]*f(x[i]); } printf("由高斯-勒让德积分公式计算得I=%lf\n",s); return 0; } 4.运行结果： 5.参考文献: [1] 谭浩强. C语言程序设计[M]. 北京：清华大学出版社，2005. [2] 秦新强. 数值逼近, 西安,2010.
数值逼近实习 题目二重积分的复化梯形公式 专业信息与计算科学 班级计算092 学号3090811072 学生薛藏朋 指导教师秦新强 2011 年
一、实验目的 1.利用Gauss-Legendre 公式计算积分 2.比较计算误差与实际误差 二、数学模型 ??? ????∈--=++≈?∑-=),(),(12)(][)]()(2)([2)(''211b a f h a b f R b f x f a f h dx x f n b a n k k ηη 三、算法 Step 1:输入等分数n Step2:输入积分上下限； Step3: 求出步长及对应个点; Step4: 由??? ????∈--=++≈?∑-=),(),(12)(][)]()(2)([2)(''211b a f h a b f R b f x f a f h dx x f n b a n k k ηη计算积分结果 Step5:将积分结果输出； 四、程序 #include
double f(double x) {double z; z=cos(x); return(z); } int main() { int i, n,m; double X[N],A[N],F=0; cout<<"请输入代数精度n"<
高斯积分点表 N = 2 节点权重 -0.5773502691896250 1.0000000000000000 0.5773502691896250 1.0000000000000000 N = 3 节点权重 -0.7745966692414830 0.5555555555555550 0.0000000000000000 0.8888888888888880 0.7745966692414830 0.5555555555555550 N = 4 节点权重 -0.8611363115940520 0.3478548451374530 -0.3399810435848560 0.6521451548625460 0.3399810435848560 0.6521451548625460 0.8611363115940520 0.3478548451374530 N = 5 节点权重 -0.9061798459386640 0.2369268850561890 -0.5384693101056830 0.4786286704993660 0.0000000000000000 0.5688888888888880 0.5384693101056830 0.4786286704993660 0.9061798459386640 0.2369268850561890 N = 6 节点权重 -0.9324695142031520 0.1713244923791700 -0.6612093864662640 0.3607615730481380 -0.2386191860831960 0.4679139345726910 0.2386191860831960 0.4679139345726910 0.6612093864662640 0.3607615730481380 0.9324695142031520 0.1713244923791700 N = 7 节点权重 -0.9491079123427580 0.1294849661688690 -0.7415311855993940 0.2797053914892760 -0.4058451513773970 0.3818300505051180 0.0000000000000000 0.4179591836734690 0.4058451513773970 0.3818300505051180 0.7415311855993940 0.2797053914892760 0.9491079123427580 0.1294849661688690
§4-4 高斯积分法及其应用 ● 由§4-3知，在计算空间等参数单元的载荷列阵及刚度矩阵时，需用到如下形式的定 积分： ηξηξd d f ?? --111 1 ),(； ζηξζηξd d d f ? ?? ---1 1111 1 ),,( 其中被积分函数f(ξ，η，ζ)一般是很复杂的，即使能够得出它的显式，其积分也是很繁的。因此，一般用数值积分来代替函数的定积分。 ● 数值积分：在积分区域内按一定规则选出一些点，称为积分点，算出被积函数f(ξ， η，ζ)在这些积分点处的值，然后再乘以相应的加权系数并求和，作为近似的积分值。 ● 数值积分的方法有多种，其中高斯积分法可以用相同的积分点数达到较高的精度， 或者说用较少的积分数达到同样的精度。 一、高斯积分法 1．一维积分的高斯公式 ● 一维积分的高斯公式 ∑? =-=n i i i f H d f 1 1 1 )()(ξξξ (4-47) 其中f(ξi )是被积函数在积分点ξi 处的数值，H i 为加数系数，n 为积分点数目。 ● 可以证明， ? 对于n 个积分点，只要选取适当的加数系数及积分点位置，能够使(4-47)式在被积分函数为不超过(2n-1)次多项式时精确成立。 ? 由于多数函数可表示成多项式形式，这种积分适应于大多数函数。 ● 例如， ? n=1时 )()(111 1 ξξξf H d f I ==?- (a) 不论f(ξ)的次数是0还是1，只需取H １=2，ξ１＝０，上式均是精确成立的。因为 ξξ10)(C C f += (b)
1 01 ()22(0)I f d C f ξξ-===?? (c) ? 当n=2时，能保证(4-47)式精确成立所允许的多项式的最高次数是3，此时，f(ξ)的通式为 332210)(ξξξξC C C C f +++= (d) 其精确积分为 20113 2 2)(C C d f I +==?-ξξ (e) 数值积分为 ) ()() ()()(3 2322210231321210122112 1 ξξξξξξξξξC C C C H C C C C H f H f H f H I i i i +++++++=+==∑= (f) 为了在C 0～C 3取任意值(包括取零值在内)时公式(f)是精确的，显然应有 221=+H H ， 02211=+ξξH H 3 2222211= +ξξH H ， 03 22311=+ξξH H 所以，应取 2,269,350,577.03 121-=- =-=ξξ 0,000,000,000.121==H H ? 同样，对于不超过五次的多项式，只要取 n=3 130.577,350,269,2ξξ=-==- 0,000,000,000.02=ξ 6,555,555,555.09 5 21===H H 9,888,888,888.09 8 3== H
⑴ 应用概率积分法进行地表任意点A(x ，y)下沉值W(x ，y)的数学表达式如下： W （X ，Y ）=W cm ·C X ·C y = m ·q · C 0 S α·C X ·C y C X =π 1[?-r x e πλ02d λ—?--r l x e πλ02 d λ] C y =π1[?-1 20r y e πλd λ—?-- 220 r L y e πλd λ] m —— 采出煤层厚度； q —— 地表下沉系数； α —— 煤层倾角； C X 、 C Y —— 分别为待求点在走向和倾向主断面上投影点处的下沉分布系数； l 、 L 一一 分别为采区拐点平移后走向长度及倾斜方向在地表的计算开采宽度； r ，r 1，r 2 —— 分别为走向、下山、上山的主要影响半经； (2) 地表任意点倾斜值i(x ，y)φ (沿指定的φ方向) i φ=i x C y cos φ+i y C X s i n φ i φ+90= —i x C y s i n φ—i y C X cos φ i x 、i y ——分别为待求点沿走向和倾向主断面上投影点处迭加后的倾斜变形值，mm/m ； (3) 地表任意点的曲率变形K φ K φ= Kx C y cos 2φ+Ky C X s i n 2φ+W i i cm y X sin2φ K φ+90= Kx C y sin 2φ+Ky C X cos 2φ-W i i Cm y X sin2φ Kx 、ky ——分别为待求点沿走向及倾向在主断面投影处迭加后的曲率值。 （4）地表任意点沿指定方向φ的水平移动值u φ u φ=u x C y cos φ+u y C X s i n φ u φ+90= —u x C y cos φ+u y C X sin φ u x 、u y ———分别为待求点沿走向和倾向在主断面投影点处迭加后的水平移动值， mm 。 (5) 地表任意点的水平变形εφ (沿指定方向) εφ=εx C y cos 2φ+εy C X s i n 2φ+ [(u x i y + u y i x )/W cm ] sin φ·cos φ εφ+90 =εx C y sin 2φ+εy C X cos 2φ- [(u x i y + u y i x )/W cm ] sin φ·cos φ εx 、εy ———分别为待求点沿走向和倾向在主断面投影点处迭加后的水平变形值， mm/m 。 上述公式即为概率积分法计算各种移动与变形的基本公式。
附錄 高斯積分法(Gaussian Quadratures) 一般積分最常見的方式是對於一個連續的函數作切割，使連續函數成離散形式作矩型面積積分，大都是取均等區間為切割的依據，而高斯積分能夠在特定的座標軸(abscissas)取值，並且給予一個加權權數(weight)能夠以少數的取點數為基礎，快速的計算出結果。其準確性必須依積分曲線的平滑程度而定，並非取點越多越精準。 簡式如下： 1()()()b N j j j a W x f x dx w f x =≈∑∫ j w :座標軸(abscissas) j x ：權數(weight) 簡單來說高斯積分法是對一個函數做定積分(definite integral)，通常是在積分定義的值域內，以正交的觀念找一些特定的點，代入多項式中，並給予相對應的權重後，相加得到逼近原始函數的積分值，其數學推導則不在本研究中詳敘。 因為依照不同的微分方程可推導出一些不同的高斯積分，其區別可依定積分之上下限區分，factor copula 主要是對常態分配或student’s t 分配做積分值域介於~?∞∞的對稱分配，可採用Gauss-Hermite 其來自於Hermite polynomials。 2 221()()()()()k n x x x k k n k f x dx e e f x dx w x e f x R x ∞∞ ?=?∞?∞????==+????∑∫∫ 12212!()[()] n k n k n w x n H x π??=
(2)()()2(2)!n n n n R x f x n = ()n H x ：Hermite 多項式