基于Romberg算法的数值积分的原理与实现
高等数值分析
数值积分方法小述一、背景数值积分方法发展的前提是在17世纪以牛顿和莱布尼茨为首的一批数学家发展起来的微积分。
在最初的研究中，求解积分的方法便是找到求解原函数的方法，得到原函数，以此为基础解决其他问题。
但是在深入的研究中，逐渐发现一些函数的原函数求解极其困难，甚至无法表示出来，是超越函数，还有的根本没有原函数，比如对于延拓函数：sin ,0()1,0xx f x xx ⎧≠⎪=⎨⎪=⎩无法求出它的原函数，这时要求它的积分就无法使用牛顿-莱布尼茨公式了，解决积分的问题便受到阻碍。
这种情况下就需要寻求一种新的求积分的方法来解决这些问题了。
数值积分方法便在数学家们的需求下发展起来。
二、发展历程等距节点的多项式插值求积法的观点最早是1676年出现在Newton 给Leibniz 的一封信中。
1711年，Cotes 在总结了牛顿的观点后，系统归纳了小于10个节点的插值求积方法，并发表了一篇相关论文。
1743年，Simpson 发表他所研究的求积方法。
但是从历史上看，对于辛普森的方法，数学家Cavalieri 和Gregory 似乎研究的更早，而且Cotes 也早就得到了这种方法。
1814年，数学王子Gauss 在研究这个问题时，通过优化那些求积节点得到一种更高精度的数值求积分方法，随后便发表了他的第一篇关于数值求积分的论文。
时间过了100多年，数学家Fejer 于1933年，将Chebyshev 点作为节点应用于数值求积分中，得到了一种新的方法。
1960年，数学家Clenshaw 和Curtis 研究得到一种更为高效的数值求积公式。
Kronrod 在1964年发表了他自己的数值求积方法，4年后的Patterson 对这种方法进行了推广，得到的方法也为世人所知。
值得一提的是Richardson 在1927年发现的外推法，当时并没有用来做数值积分问题。
而数学家Romberg 在1955年将它应用到数值积分上，取得不小的成果。
数值分析MATLAB编程——数值积分法
数值分析MATLAB编程——数值积分法1、调用函数--f.Mfunction y=f(x)%------------------------------------------------------------函数1 y=sqrt(4-sin(x)*sin(x));%------------------------------------------------------------函数2 %y=sin(x)/x;%if x==0% y=0;%end%------------------------------------------------------------函数3 %y=exp(x)/(4+x*x);%------------------------------------------------------------函数4 %y=(log(1+x))/(1+x*x);2、复合梯形公式--tixing.M%复合梯形公式clear alla=input('请输入积分下限:');b=input('请输入积分上限:');n=input('区间n等分：');h=(b-a)/n;x=a:h:b;T=0;for k=1:n;T=0.5*h*(f(x(k))+f(x(k+1)))+T;endT=vpa(T,8)3、复合Simpson公式--simpson.M%复合Simpson公式clear alla=input('请输入积分下限:');b=input('请输入积分上限:');n=input('区间n等分：');h=(b-a)/n;x=a:h:b;S=0;for k=1:n;xx=(x(k)+x(k+1))/2;S=(1/6)*h*(f(x(k))+4*f(xx)+f(x(k+1)))+S;endS=vpa(S,8)4、Romberg算法--romberg.M%Romberg算法clear alla=input('请输入积分下限:');b=input('请输入积分上限:');n=input('区间n等分：');num=0:n;R=[num'];h=b-a;T=h*(f(a)+f(b))/2;t(1)=T;for i=2:n+1;u=h/2;H=0;x=a+u;while x<b;H=H+f(x);x=x+h;endt(i)=(T+h*H)/2;T=t(i);h=u;endR=[R,t'];for i=2:n+1for j=n+1:-1:1if j>=it(j)=(4^(i-1)*t(j)-t(j-1))/(4^(i-1)-1);elset(j)=0;endendR=[R,t'];endR=vpa(R,8)R(n,n)5、变步长算法（以复化梯形公式为例）--tixing2.M%复合梯形公式，确定最佳步长format longclear alla=input('请输入积分下限:');b=input('请输入积分上限:');eps=input('请输入误差：');k=1;T1=(b-a)*(f(a)+f(b))/2;T2=(T1+(b-a)*(f((a+b)/2)))/2; while abs((T1-T2)/3)>=epsM=0;n=2^k;h=(b-a)/n;T1=T2;x=a:h:b;for i=1:n;xx=(x(i)+x(i+1))/2;M=M+f(xx);endT2=(T1+h*M)/2;k=k+1;endT=vpa(T2,8)n=2^k。
【推荐】数值计算方法：第5章-数值微分与数值积分.ppt
20
（1）插值型求积公式
2.由下列列表函数求L-插值多项式
x0
x1 --- xi-1
xi
xi+1
---
xn
f(x0) f(x1) --- f(xi-1) f(xi) f(xi+1) --- f(xn)
21
称为插值型求积公式，
称为求积节点， 称为求积系数，其和
22
求积系数 通过插值基函数
这称为梯形公式；
a
b
图1 梯形公式
几何意义：用梯形面积 代替f(x)作为曲边的曲边 梯形面积。
25
这称为Simpsion公式。
a
b
图2 Simpson公式
几何意义：用抛物线 作曲边的曲边
梯形面积代替f(x)作 为曲边的曲边梯形面26积。
这称为Cotes公式。
求积公式的误差（余项）
27
28
例 5. 1 分别用梯形公式、Simpson公式计算定积分
时的步长h/2就是合适的步长
6
例：
f(x)=exp(x)
h
f’(1.15) R(x)
h
f’(1.15) R(x)
0.10 0.09 0.08 0.07 0.06
3.1630 -0.0048
3.1622 3.1613 3.1607 3.1600
-0.0040 -0.0031 -0.0025 -0.0018
17
构造数值积分公式的基本思想： 由积分中值定理知，在积分区间
成立
内存在一点ξ，
问题：点ξ的具体位置一般是不知道的，因而难以 准确算出 的值，怎么办？
只要对平均高度 一种数值求积方法.
提供一种算法，相应地便可获得
第2章数值微分和数值积分
f '( x) D(h) O ( h) 2 f '( x) D(h / 2) O(h / 2)
f '( x) D(h) 2 f '( x) 2D(h / 2) f '( x) D(h / 2) D(h) D(h / 2)
例：
f(x)=exp(x)
h 0.10 0.09 f’(1.15) 3.1630 3.1622 R(x) -0.0048 -0.0040 h 0.05 0.04 f’(1.15) R(x) 3.1590 3.1588 -0.0008 -0.0006
f
误差
(k )
( x) Ln ( x)
(k )
f ( n1) ( ) Rn ( x) n ( x) f ( x) Ln ( x) (n 1)! k ( n 1) d f ( ) (k ) Rn ( x) k n ( x) dx (n 1)!
1 h2 f '( x0 ) L '2 ( x0 ) 3 f ( x0 ) 4 f ( x1 ) f ( x2 ) f '''( ) 2h 3 1 h2 f '( x1 ) L '2 ( x1 ) f ( x0 ) f ( x2 ) f '''( ) 2h 6 1 h2 f '( x2 ) L '2 ( x2 ) f ( x0 ) 4 f ( x1 ) 3 f ( x2 ) f '''( ) 2h 3 1 h2 (4) f ''( x0 ) L ''2 ( x0 ) 2 f ( x0 ) 2 f ( x1 ) f ( x2 ) [ hf '''(1) f (2 )] h 6 1 h2 (4) f ''( x1 ) L ''1 ( x2 ) 2 f ( x0 ) 2 f ( x1 ) f ( x2 ) f ( ) h 12 1 h2 (4) f ''( x2 ) L ''2 ( x2 ) 2 f ( x0 ) 2 f ( x1 ) f ( x2 ) [hf '''(1 ) f (2 )] h 6 2 Taylor展开分析，可以知道，它们都是 O(h )
gauss-kronrod 自适应数值积分算法
一、介绍gauss-kronrod 自适应数值积分算法是一种用于解决数值积分问题的高精度算法。
它基于高斯-克罗德积分公式，能够在较少采样点的情况下取得较高的积分精度。
该算法被广泛应用于科学计算、工程领域和数学建模中，能够有效地解决复杂函数的数值积分问题。
二、算法原理1. 高斯-克罗德积分公式高斯-克罗德积分公式是一种基于重要节点的数值积分方法。
通过在特定节点上取样，可以将积分问题转化为对这些节点处函数值的加权求和，从而得到积分的近似解。
高斯-克罗德积分公式具有较高的精度和稳定性，在一定程度上可以避免数值积分中的波动和误差。
2. 自适应策略gauss-kronrod 算法在高斯-克罗德积分公式的基础上，加入了自适应策略。
它能够根据积分函数的特性和精度要求，动态调整采样点的位置和数量，使得在给定精度下能够取得较高的积分精度。
这种自适应策略在面对复杂的积分函数时尤为有效，能够减少不必要的计算量，提高积分的效率和准确性。
三、算法实现1. 积分区间划分在使用 gauss-kronrod 算法进行数值积分时，首先需要对积分区间进行适当的划分。
与传统的数值积分算法不同的是，gauss-kronrod 算法能够根据积分函数的特性，自适应地调整划分的方式和密度，从而更有效地逼近积分的精确值。
2. 采样点选择gauss-kronrod 算法基于高斯-克罗德积分公式，采用了一组经过优化的重要采样点。
这些采样点的位置和权重经过精心设计，能够在保证积分精度的前提下，减少采样点的数量，提高计算效率。
3. 自适应策略调整在进行数值积分计算时，gauss-kronrod 算法会根据当前精度的要求和实际的积分情况，动态调整采样点的位置和数量。
通过在积分过程中不断进行自适应调整，能够有效地减少不必要的计算量，提高计算效率，并且保证积分的精度。
四、算法优势1. 高精度gauss-kronrod 算法基于高斯-克罗德积分公式，能够在少量采样点的情况下取得较高的积分精度。
用MATLAB计算某些区域上的二重积分
用ＭＡＴＬＡＢ计算某些区域上的二重积分摘要：本文研究某些区域上二重积分数值积分公式的构造及用数学软件MATLAB 实现所构造的数值积分的计算，通过MATLAB 用所得公式计算某些典型的二重积分。
主要工作包括：将定积分数值计算的几个公式推广到二重积分，编制MATLAB 程序，最后通过具体的数值算例进行精度比较，从中选出精度高的二重积分的计算公式，并用公式计算一些典型的二重积分。
关键词：二重积分；数值计算；插值多项式；求积公式 1.引言二重积分的计算在科学计算中起着重要的作用，关于矩形区域上的二重积分的计算一般都是化重积分为累次积分，然后借助定积分已有的数值积分计算公式推导出，MATLAB 已经有这些计算公式的相应的函数，但是往往我们建模得到的二重积分的积分区域都不是矩形区域，对于一般的非矩形区域的二重积分，直接用MATLAB 是无法计算的。
又当被积函数比较复杂，无法用初等函数表示或求其原函数很困难时，就只能求积分的数值解。
若),(y x f 在D 上连续，二重积分y x y x f Dd d ),(⎰⎰存在且为一确定的常数，这个数值与),(y x f 的结构、D 的几何形状有关，二重积分计算的基本途径是在一定条件下化为二次积分，本文研究的某些区域的二重积分，要求二重积分在该区域上能化为二次积分。
二重积分的存在性[1]：),(y x f 在闭区域D 上连续，则(,)Df x y dxdy ⎰⎰必存在。
定理[1]：若),(y x f 在闭区域D {;()()}≤≤≤≤a x b c x y d x 上连续，且()c x 、()d x 在[,]a b 上连续，则：()()(,)(,)b d x ac x DI f x y dxdy f x y dxdy ==⎰⎰⎰⎰上式右端是一个先对x 后对y 的二次积分：先把),(y x f 看作x 的函数，在区间[(),()]c x d x 上对x 计算定积分（这时y 看作常数），把得到的结果（是y 的函数）再在[,]a b 上对y 计算定积分即为二重积分。
采用数值重积分法精确求解辊底炉内角系数
冶金能源
Sept. 2019
ENERGY FOR METALLURGICAL INDUSTRY
31
采用数值重积分法精确求解棍底炉内角系数
蒋彦军
（宝钢工程技术集团有限公司）
摘要基于Romberg积分公式，使用计算机递归算法，建立直接求解重积分的数值方法, 精确求解辐底炉内各辐射面之间的角系数。该方法精度可控、积分速度快，可以直接推广至 其他炉型，为三维辐射换热的精确求解奠定基础。 关键词角系数数值重积分辐底炉辐射换热 文献标识码：A 文章编号：1001 -1617 （2019） 05 -0031 -06
角系数的确定方法主要有：（1）计算公式 及图表法；（2）图解法；（3）实验测定法；（4） 数值计算法。数值计算是求解该问题可行的方
收稿日期=2019 - 05 - 05 蒋彦军（1983 -）,工程师；201900上海市宝山区。
法，主要包括有限差分法、有限元法⑴、边界 元法⑵、蒙特卡洛法⑶和矢量法等⑷。这些方 法都有一定的不足，如有限元法计算精度取决于 离散网格数量，不易控制计算精度；蒙特卡洛法 为了达到统计模拟的目的，要求每个单元面积发 射的能束数量要有足够的数目，且存在很多重复 计算，需要耗费大量的计算时间⑸。矢量法一 般情况下只适用于简单结构之间的角系数求 解叫
（衍，光2，・"_1）和仇（衍，光2,…省一1 ）的解析。 两个辐射面相交或接近的情况下会出现奇异
积分，数值求解时部分积分出现无法收敛、收敛
慢等情况，可将相交处的面域适当微量剪除，剪
除部分可按1取角系数值，计算直接交换面积后
叠加到积分结果里。
部分辐射面远大于相对间距或窄条型面域会
出现收敛慢的情况，合理重新排列各维度的先后
若干数值积分的计算方法
1 若干数值积分的计算方法 黄海琼 （广西民族大学数计学院04数本1班 南宁 530006）
摘 要: 本文讨论了若干数值积分的计算方法。在一维情形下，介绍了Newton-Cotes公式，Gauss型等求积法则； 在二维情形下, 主要介绍了二元Newton-Cotes积分方法。最后，对几类数值积分方法及其数值实验进行比较评述。 关键词: 牛顿-柯特斯公式；Gauss型求积法则；二元数值积分；数值实验
Some Computational Methods of numerical integration Huang Haiqiong (College of Mathematics and Computer Science,Guangxi University for Nationalities, Nanning 530006) Abstract: In this paper, some computational methods about numerical integration are discussed. under the univariate situation, the quadrature rule of Newton-Cotes formula, Gauss formula and so on is introduced. Under the two-dimensional situation, it mainly introduced the dual Newton-Cotes integral method. Finally, the numerical integration methods and numerical experiment were discussed. Key word: Newton-Cotes formula; Gauss integration principle; dual numerical integration; numerical experiment.
数值积分实用PPT课件PPT课件
微积分学中，积分计算是利用 Newton – Leibniz
公式：
来计算的。
例，某气体由温度 T1 加热到 T2 时所需热量 Q 可由下式表示：
Q
T1 T2
Cp
.mdT
Cp .m 该气体的摩尔定压热容
不知道该气体的 与CTp的.m函数关系式，而实验测得该气体的
系数据如下表所示。
C p .m
与 T 的关
T/℃
25 100 150 200 250 300 350 400 450 500 600 700 800
32 k0
f
( x ) 12
k1
k0
4
f
( xk1 ) 2
n1
n1
32 f ( x ) 14 f ( x )
k0
k3 4
k 1
k
第27页/共39页
龙贝格算法
第28页/共39页
数值方法中常利用一序列{ F1、F2、…、Fk、…} 去逼近精确值，然后在理论上给出序列F的误差估计。
新思路：
能否在某种理论（截断误差估计）基础上，通过简单方法，在序列
值，而这些值又是成倍增加的，所以计算工作量较大。
第37页/共39页
程序例子：P207
作 业：
将P207的程序改为求解积分
eps=0.000001 结果（0.1115718）
1x
0 4 x2 dx
第38页/共39页
感谢您的观看。
复化求积分
)
实际计算中的递推公式为 b−a T1 = [ f (a ) + f (b )] 2 1 b − a n −1 b−a T2 n = Tn + ∑ f (a + (2 j + 1) 2n ), 2 2n j = 0
n = 1, 2, ⋯
直到 | T2n − Tn |≤ ε 为止，T2n作为积分的近似值。
假定f ''( x )在[a , b]上变化不大, 即有f ''(η1 ) ≈ f ''(η 2 )， 于是得
I − Tn ≈4 I − T2 n
1 1 ∴ I ≈ T2 n + (T2 n − Tn ) 或 I − T2 n ≈ (T2 n − Tn ) 3 3 1 I 当 T2 n − Tn ≤ ε 时 ， − T2 n ≤ ε ≤ ε 。 3 数值分析
数值分析
数值分析
同理可得变步长复化柯特斯公式
实际计算过程如下： 实际计算过程如下： T1 → T2 ↓ S1 → T4 ↓ S2 ↓
C1 →
→ →
T8 ↓ S4 ↓
C2 ↓ R1
→ →
T2 n =
1 1 Tn + H n 2 2
h2 h4
1 S n = ( 4T2 n − Tn ) 3
42 S 2n − S n Cn = 42 − 1 43 C 2n − C n Rn = 43 − 1
1 2 已 知 S n = Tn + H n 3 3 1 1 又 有 T2 n = Tn + H n 2 2 两式联立解得： 1 2 1 S n = Tn + （2T2 n − Tn） （4T2 n − Tn） = 3 3 3
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
基于Romberg算法的数值积分的原理与实现
杨俊峰’ 李静
(1．贵州大学电信学院 贵州贵阳 550003； 2．菏泽学院计算机系 山东菏泽 27401 5)
摘要：本文首先简单地介绍了数值积分的基本思想，然后阐述了复合梯形公式、Richardson
外推算法及Romberg算法，最后给出
了实现Romberg算法的通用函数。
关键词：数值积分 复合梯形公式 Richardson外推算法 Romberg算法
中图分类号：G424 文献标识码：A 文章编号：1673 0534(2o07)03(b)一0059—02
在实际的计算中，常常需要计算定积分
l f(x)dx的值，若知道f(x)的原函数V(x)，
则由牛顿 莱布尼兹公式就把求定积分的问
题转化成了求函数F(b)～F(a)的值了。但是，
实际上很多被积函数的原函数不能用初等函
数表示出来，有的即使可用初等函数表出但
表达式较复杂，再者在实际应用中，有时被积
函数是以表格或图形的形式给出，因此对千
这些问题，需要建立数值积分的方法才能加
以解决。
1数值积分的基本思想川
定积分的几何意义是曲边梯形的面积，
而在高等数学上，定积分的定义是求和式的
极限，和式中每一项则是小曲边梯形面积的
近似值 小矩形的面积。计算中并不能取
极限这个无限过程，这就需用有限个小矩形面
积之和来近似代替曲边梯形面积，这实际上
是用阶梯函数来近似代替被积函数。例如构
造一个简单函数(如多项式) (x)来近似代替
一．
被积分函数，(X)，然后通过求1 ( )出
求得I八 )dx的近似值。为了精确计算，往
往把积分区间细分，然后再在每个小区间内
用这个简单函数代替被积函数并积分以求得
积分的近似值，这就是数值积分的基本思想。
2 Romberg方法数值积分的相关原理[2】【。1
2 1复合梯形公式
由于梯形公式
}、
l八 )ch 【 ((1 ．，(^1】，误差较大，一
般无法使用。下面讨论改进的梯形公式即复
合梯形公式，在这里为了说明问题的方便，将
积分区间【a，bln等分，节点
h 一口
：a+ih，h： ，i=1,2，3⋯，，rl，对每
一个小区间Ix．，Xi+l】采用梯形公式，就得到如
下的复合梯形公式：
m)=ff(x Σn-1 )+ 川 】：
．，’(口)+2艺厂(口+ )+厂(6)J，h表示把
积分区间Ia，b】n等分，这样T、h )就可表示把
积分区间Ia，bl2n等分，则可 推导出
九 h h
x ){
汶个关系给}}j了 积分 问的分点增加后的
计算公式与原计算公式之间的关系，即当节
h
点增加一倍计算复合梯形公式T(三)时，只需
把原来复合梯形公式T(h)的值除以2，再加上
’， ， 、h
新节点上的函数值 + 之和与新步长
的乘积。这样既可使计算结果更精确而且节
约计算量，也可由此来判断是否达到要求的
精度。
2．2 Richardson外推算法
在高等数学中，一般只考虑收敛性，并
不考虑收敛速度问题，但在实际的数值计算
中，收敛的快慢显得特别重要，因此要研究加
速收敛的技巧。
2．2．I加速收敛技巧的基本思想
在实际的数值计算中，常常用一系列数
{Fk}去逼近准确值F ，若固定k，则F 固
定，随之误差F‘F 也就固定。因此要提高
计算的精确度，就必须计算F +．，F + ，
F ，⋯，现在需要的是不去计算新值F ，
F +，，F + ， ⋯，而是由已知信息F．，F，，
F ，⋯ ，F ．，F 的线性组合，得出新的
序列{F 1 ，而使F1．F+的误差较小， 亦即
f J }比{F， }更快地收敛于F ，这是加速收敛技
巧的基本思想。
2．2．2 Richardson外推算法
假设有一命题F ，用一个步长为h的函
数F(h)去逼近(这里F 与h无关)，则Richar
dson外推算法由如下的递推关系来实现
0)F。( )= h)
(2) I(h)= ，
m：1，2，3⋯ ，P 是与h无关的常数，q为满
足1一q ≠O(m：1．2．3⋯ ．)的适当正数。
Richardson外推算法流程图如下
(2)F IlqI、严=。。 )，=—=◆： ‘S)F qh) ⋯” ’
(4)F{(q—h)／一 ⋯ ⋯
2，3 Romberg算法
该算法是在计算梯形和序列的基础上应
用了线性外推的加速方法， 由此构成的一种
具有超线性收敛的数值积分法， 即把
2 1
，9 代入Richardson#l-推算法的递
推公式得到一个更快地收敛于I )f缸的新
序列，记为 ，称为Romberg序列，其
中k表示将积分区间2 等分，m表示外推m
1次，即m=I时表示不外推，如T．(。 表示在
【a，bl上的梯形公式，T，( 表示把区间2 等
分的复合梯形公式，而T，‘ 表示对TI(U~T．(
外推所得的公式。这样Romberg算法的具体
描述如下：
(1)先求出按梯形公式所得的积分值
r ，=—b
—：
-
：一a[，((f)+，(6)】
(2)把区间2等分，求出两个小梯形的面积
之和，记为T ‘”，即
= 『，(小 2，{ 1
再由外推法可得T =
(3)把区间再等分(即2 等分)，得复合梯形
公式T．‘ ，由T． ’与T1‘ 外推可得
T ”=
，
T =兰
如此若己计算出2 等分的复合梯形公式
T，‘ ，则由Richardson外推算法就可构造出
新序列
T‘ ，： ， m = lf2，．．．．，．p，
m+{ 一 1 ’ ¨1 ⋯ ’ P’
k：1,2一 ，P—m+l，最后求得T (O 。
(4)当T T ： 或lT ! 一T l≤ ，(
为给定的任意小的正数)就停止计算，否则返
N(3)，计算T1‘⋯ 。
根据匕算法可以得到如下的具体流程：
(1) TI'。
(2)
(4)
。 r。◆
f51 T： ⋯ ⋯
3 Romberg算法的实现程序
根据上面对Romberg算法原理的分析与
描述，不难得到只要把T ( ’用二维数组来表
示基本上就完成了该算法到程序的实现。下
面给出，只要知道被积函数及其下限a 上限
b和要达到的精度 ，都可以计算数值积分的
通用函数(其实该函数的返回值就是要求的数
值积分值)，具体函数如下
double Romberg(double(*f)(double)，
double a，double b，double )
{int P，i，k，m；
double TINIIN】，x；
T[1l[0l=((b～a)／2) (( f)(a)+(+f)(b))}／／
初始化TI1]10J
p=1；
while(I)
{x=0；__for(i~1；i<=pow(2，P 1)；i++)／／计算
T【l】【p】
x=x+f(a+(2+i 1)+((b a)／pow(2，p)))；
x=(x}((b a)／pow(2，P 1))+T【l】【p—
l1)／2；
T【l】【p】=x；
k=p；
for(m=l；m<=P；m++) ／／利用
Richardson外推法计算T【m+l】【k l】
{T【m+l】【k 1]=(pow(4，m)}T【m】【k卜T
【m】【k 1])／(pow(4，m)1)；
k ；}
if(fabs(T[p+l】【O卜T【p】【O】)<)
{return(T[p+lⅡO】)；／／返回满足题目要
求精度的计算结果T【p+l】【O】
breakI}
else
p++ ；
}}
4结语
用Romberg算法计算数值积分时，一般
计算到p=4或5即能达到要求，该算法简单
且收敛速度快，也节省存储容量，但是对函
数的光滑性要求较高。
参考文献
⋯1 王能超．数值分析简明教程．高等教育出版
60 科技咨询导报Science and Technology Consulting Herald
社，2002．8．
【2】何伟保，张民选．数值分析．贵州科技出版
社，2003．8．
【3】沈剑华．数值计算基础．同济大学出版社，