龙贝格积分
龙贝格求 积分
龙贝格(Romberg ）求积法1.算法理论Romberg 求积方法是以复化梯形公式为基础，应用Richardson 外推法导出的数值求积方法。
由复化梯形公式 )]()(2)([2222b f h a f a f h T +++=可以化为)]()]()([2[212112h a f h b f a f hT +++=＝)]([21211h a f h T ++一般地，把区间［a,b ］逐次分半k －1次，（k ＝1，2,……,n)区间长度（步长)为kk m a b h -=，其中mk ＝2k －1。
记k T ＝)1(k T由)1(k T ＝]))12(([21211)1(1∑=---++km j k k k h j a f h T 从而⎰badxx f )(＝)1(kT－)(''122k f h a b ξ- (1)按Richardson 外推思想，可将（1）看成关于k h ，误差为)(2k h O 的一个近似公式，因而，复化梯形公式的误差公式为⎰badxx f )(－)1(k T ＝......4221++kkh K h K ＝∑∞=12i i k i h K （2）取1+k h ＝k h 21有 ⎰ba dx x f )(－)1(1+k T ＝∑∞=+121221i ik ii hK （3）误差为)(2jh O 的误差公式 )(j kT＝)1(-j kT＋141)1(1)1(------j j k j k T T2。
误差及收敛性分析（1）误差，对复化梯形公式误差估计时，是估计出每个子区间上的误差,然后将n 个子区间上的误差相加作为整个积分区间上的误差。
(2）收敛性,记h x i =∆，由于∑=++=ni i i n x f x f h f T 01))]()([2)(＝))()((21101∑∑-==∆+∆n i ni i i i i x x f x x f上面两个累加式都是积分和，由于)(x f 在区间],[b a 上可积可知，只要],[b a 的分划的最大子区间的长度0→λ时,也即∞→n 时,它们的极限都等于积分值)(f I 。
梯形法的递推化和龙贝格算法
梯形法的递推化和龙贝格算法
梯形法（Trapezoid Rule）是一种数值积分的方法，用于计算函数在给定区间上的定积分。
梯形法的递推化可以通过将区间等分为若干个小梯形来实现。
具体步骤如下：
1. 将给定区间[a, b]分成n个等距的子区间，每个子区间的宽度为h = (b - a) / n。
2. 计算每个子区间的矩形面积（等于底边长度乘以高度），然后将所有子区间的矩形面积相加，得到梯形法的近似积分值。
梯形法的递推公式可以表示为：
I = h/2 * [f(a) + 2*f(a+h) + 2*f(a+2h) + ... + 2*f(a+(n-1)h) + f(b)] 其中，f(x)为要积分的函数。
龙贝格算法（Romberg Integration）是一种数值积分的方法，可以通过多次应用梯形法来逐步提高积分结果的精度。
龙贝格算法的具体步骤如下：
1. 利用梯形法计算区间[a, b]上的第一次近似积分值T(1,1)。
2. 利用递推公式T(m, 1) = 1/2 * [ T(m-1, 1) + h(m-1) *
Sigma(2^(m-2) * f(a + (2k-1) * h(m-1))), k=1 to 2^(m-2)]，计算T(m, 1)，其中h(m-1)为区间[a, b]的步长。
3. 计算T(m, n) = T(m, n-1) + 1 / (4^(n-1) - 1) * ( T(m, n-1) - T(m-1, n-1) )，其中n > 1，m > n。
4. 重复步骤3，直到达到所需的精度要求。
龙贝格算法通过递归和递推来不断提高积分结果的精度，可以较快地得到比较准确的近似积分值。
龙贝格算法
龙贝格积分1. 算法原理采用复化求积公式计算时，为使截断误差不超过ε，需要估计被积函数高阶导数的最大值，从而确定把积分区间[]b a ,分成等长子区间的个数n 。
首先在整个区间[]b a ,上应用梯形公式，算出积分近似值T1；然后将[]b a ,分半，对 应用复化梯形公式算出T2；再将每个小区间分半，一般地，每次总是在前一次的基础上再将小区间分半，然后利用递推公式进行计算，直至相邻两个值之差小于允许误差为止。
实际计算中，常用ε≤-n n T T 2作为判别计算终止的条件。
若满足，则取n T f I 2][≈；否则将区间再分半进行计算，知道满足精度要求为止。
又经过推导可知，∑=-++=ni i i n n x x f h T T 112)2(221，在实际计算中，取kn 2=，则k a b h 2-=，112)1*2(2++--+=+k i i ab i a x x 。
所以，上式可以写为∑=++--+-+=+kk i k k ab i a f a b T T 211122)2)12((2211k开始计算时，取())()(21b f a f ab T +-=龙贝格算法是由递推算法得来的。
由梯形公式得出辛普森公式得出柯特斯公式最后得到龙贝格公式。
根据梯形法的误差公式，积分值n T 的截断误差大致与2h 成正比，因此步长减半后误差将减至四分之一，即有21114n n T T -≈-将上式移项整理，知2211()3n n n T T T -≈-由此可见，只要二分前后两个积分值n T 和2n T 相当接近，就可以保证计算保证结果计算结果2n T 的误差很小，这种直接用计算结果来估计误差的方法称作误差的事后估计法。
按上式，积分值2n T 的误差大致等于21()3n n T T -，如果用这个误差值作为2n T 的一种补偿，可以期望，所得的()222141333n n n n n T T T T T T =+-=-应当是更好的结果。
龙贝格积分公式
龙贝格积分公式
龙贝格积分公式，是数学中常见的一种积分方法。
它通过分割区间，将被积函数转化为$Polynomial$（多项式）的形式，并通过加权平均的方式求出积分值。
这种方法被广泛应用于科学计算领域，如物理、化学等。
龙贝格积分公式是从重复使用$Simpson$和$Mid-point$公式推导而来的。
该公式基于分治思想，将整个区间分成若干个子区间，并对每个子区间进行逐层递推，最终得出整个区间的积分值。
在推导龙贝格积分公式时，需要利用“函数逼近”的思想，即将被积函数转化为多项式的形式。
这样可以大大简化计算，减小误差，并提高计算精度。
公式的具体计算过程如下：
假设被积函数为$f(x)$，积分区间为$[a,b]$，将积分区间均分成$2^n$个小区间，在每个小区间上做$Simpson$公式近似积分，得到$S_{2^n}$，即：
$$S_{2^n}=\frac{4^nS_{2^{n-1}}-S_{2^{n-1}}}{4^n-1}$$
其中，$S_{2^n}$为$n$级逼近值，$S_{2^{n-1}}$为$n-1$级逼近值。
根据上式，可得$S_{2^1}$，然后再计算$S_{2^2}$，$S_{2^3}$，以此类推，递归地计算$n$级逼近值，直到计算所得值与精确值的差别小于预先设定的精度要求为止。
龙贝格积分公式没有强制要求$f(x)$连续可微，又由于是基于函数逼近的方式进行积分，精度高且计算速度快，因此被广泛应用。
总之，龙贝格积分公式是一种有效的求解复杂积分问题的方法，在处理高维积分时，具有更大的优势。
4.4龙贝格求积公式
4 1 T1 ( k 1) T0 ( k ) T0 ( k 1) 3 3 16 1 T2 ( k 1) T1 ( k ) T1 ( k 1) 15 15 64 1 T3 ( k 1) T2 ( k ) T2 ( k 1) 63 63
k 1, 2 ,
1 h n1 Tn f ( x 1 ) j 2 2 j 0 2 1 b a n 1 1 Tn f (a ( j )h) 2 2n j 0 2 1 b a n 1 ba Tn f ( a ( 2 j 1) ) 2 2 n j 0 2n
Romberg算法的代 数精度为m的两倍 Romberg算法的收敛 阶高达m+1的两倍
T0 (0 ) T0 (1) T0 ( 2 ) T0 ( 3) T1 (0 ) T1 (1) T1 ( 2 ) T2 (0 ) T2 (1)
T3 (0 )
龙贝格积分
例：计算
I T2n 1 I Tn 4
令
--------(7)
即 当然
Cn C2 k 1 T2 (k 1) C 2 n T2 (k )
--------(8)
同样由复合Cotes公式的余项
I C2 n 1 (C2 n Cn ) 63
64 1 64 1 得 I C2 n Cn T2 ( k ) T2 ( k 1) 63 63 63 63
43 C 2n C n Rn 3 4 1
Romberg
T1 = T0( 0 )
T2 = T0( 1 )
<?
算法： T4 = T0( 2 )
T8 = T0
(3)
S1 = T1( 0 ) S2 = T1
龙贝格算法求数值积分程序说明及流程图
In(k,2)=4*In(k+1,1)/3-In(k,1)/3;
end
flag(2)=i+3;
for k=flag(3)+1:i+2
In(k,3)=16*In(k+1,2)/15-In(k,2)/15;
end
flag(3)=i+2;
for k=flag(4)+1:i+1
end
In(i+5,1)=In(i+5,1)/2^(i+5)+In(i+4,1)/2;
In(i+4,2)=4*In(i+5,1)/3-In(i+4,1)/3;
In(i+3,3)=16*In(i+4,2)/15-In(i+3,2)/15;
In(i+2,4)=64*In(i+3,3)/63-In(i+2,3)/63;
陆韶琦 3110000441
程序说明：本程序计算
数值积分值。
b sinx
a x
dx，程序初始要求输出需要得到的精度，最后输出得到
输入精度 eps，积分上下限 b,a
流程图：
定义函数 f(x)=
1
x=0;
sin(x)/x 其他。
计算 Ti(0<i<5)
T2^(i+1)=T2^i/2+
2j−1 b−a
2^i
j=1 f(
2 i +1
+ )
计算 Si（0<i<4）
Si=4Ti+1/3-Ti/3
计算 Ci(0<i<3)
龙贝格算法
龙贝格算法一、问题分析1、1龙贝格积分题目要求学生运用龙贝格算法解决实际问题（塑料雨篷曲线满足函数y（x)=l sin （tx），则给定雨篷得长度后，求所需要平板材料得长度）.二、方法原理2、１龙贝格积分原理龙贝格算法就是由递推算法得来得。
由梯形公式得出辛普生公式得出柯特斯公式最后得到龙贝格公式.在变步长得过程中探讨梯形法得计算规律.设将求积区间［ａ，b]分为ｎ个等分，则一共有n+1个等分点，n.这里用表示复化梯形法求得得积分值，其下标n 表示等分数。
先考察下一个字段[]，其中点,在该子段上二分前后两个积分值显然有下列关系将这一关系式关于k从０到n－１累加求与，即可导出下列递推公式需要强调指出得就是,上式中得代表二分前得步长，而梯形法得算法简单，但精度低,收敛速度缓慢,如何提高收敛速度以节省计算量，自然式人们极为关心得.根据梯形法得误差公式，积分值得截断误差大致与成正比，因此步长减半后误差将减至四分之一,既有将上式移项整理,知由此可见，只要二分前后两个积分值与相当接近，就可以保证计算保证结果计算结果得误差很小，这种直接用计算结果来估计误差得方法称作误差得事后估计法。
ﻩ按上式,积分值得误差大致等于，如果用这个误差值作为得一种补偿，可以期望,所得得应当就是更好得结果。
ﻩ按上式，组合得到得近似值直接验证，用梯形二分前后得两个积分值与按式组合，结果得到辛普生法得积分值。
再考察辛普生法。
其截断误差与成正比.因此,若将步长折半,则误差相应得减至十六分之一。
既有由此得不难验证,上式右端得值其实就等于,就就是说,用辛普生法二分前后得两个积分值与，在按上式再做线性组合,结果得到柯特斯法得积分值，既有重复同样得手续，依据斯科特法得误差公式可进一步导出龙贝格公式应当注意龙贝格公式已经不属于牛顿—柯特斯公式得范畴.在步长二分得过程中运用公式加工三次，就能将粗糙得积分值逐步加工成精度较高得龙贝格,或者说,将收敛缓慢得梯形值序列加工成熟练迅速得龙贝格值序列，这种加速方法称龙贝格算法。
龙贝格求积公式
Romberg积分公式正是由此思想产生 积分公式正是由此思想产生
类似于梯形加速公式的处理方法，得到： 类似于梯形加速公式的处理方法，得到： 加速公式的处理方法
2 1 4 S2 n − Sn I ≈ S2 n + 2 ( S2 n − Sn ) = 2 4 −1 4 −1 2 4 S2 n − Sn Simpson加速公式： C n = 加速公式 加速公式： 2 4 −1 3 1 4 C2n − Cn I ≈ C 2 n + 3 (C 2 n − C n ) = 3 4 −1 4 −1
求出f(1/4)
f(3/4) 进而求得
T4=1/2{T2+1/2[f(1/4)+f(3/4)]} =3.131176471 S2=1/3(4T4-T2)=3.141568628 C1=1/15(16S2-S1)=3.142117648 计算f(1/8) f(3/8) f(5/8) f(7/8)进而求得 T8=1/2{T4+1/4[f(1/8)+f(3/8)+f(5/8)+f(7/8)]} =3.138988495 S4=1/3(4T3-T4)=3.141592503 C2=1/15(16S4-S2)=3.141594095 R1=1/63(64C2-C1)=3.141585784
梯形值序列
4个积分值序列： 个积分值序列： 个积分值序列
T2k
S 2k =
C 2k =
4T2k +1 − T2k 4 −1 2 4 S2k +1 − S2k
4 −1
2
2k
2k 2k
R2k =
4 C 2k +1 − C 2k
3
4 −1
3
龙贝格积分算法实验
实验题目2 Romberg 积分法摘要考虑积分()()b aI f f x dx =⎰欲求其近似值，可以采用如下公式： （复化）梯形公式 11[()()]2n ii i hT f x f x-+==+∑2()12b a E h f η-''=-[,]a b η∈ （复化）辛卜生公式 11102[()4()()]6n i i i i hS f x f x f x -++==++∑4(4)()1802b a h E f η-⎛⎫=- ⎪⎝⎭ [,]a b η∈ （复化）柯特斯公式 111042[7()32()12()90n i i i i hC f x f x f x -++==+++∑31432()7()]i i f xf x +++6(6)2()()9454b a h E f η-⎛⎫=- ⎪⎝⎭[,]a b η∈ 这里，梯形公式显得算法简单，具有如下递推关系121021()22n n n i i h T T f x -+==+∑因此，很容易实现从低阶的计算结果推算出高阶的近似值，而只需要花费较少的附加函数计算。
但是，由于梯形公式收敛阶较低，收敛速度缓慢。
所以，如何提高收敛速度，自然是人们极为关心的课题。
为此，记0,k T 为将区间[,]a b 进行2k等份的复化梯形积分结果，1,k T 为将区间[,]a b 进行2k等份的复化辛卜生积分结果，2,k T 为将区间[,]a b 进行2k等份的复化柯特斯积分结果。
根据李查逊（Richardson ）外推加速方法，可得到1,11,,0,1,2,40,1,2,41m m k m km k m k T T T m -+-=-⎛⎫=⎪=-⎝⎭可以证明，如果()f x 充分光滑，则有,lim ()m k k T I f →∞= （m 固定）,0lim ()m m T I f →∞=这是一个收敛速度更快的一个数值求积公式，我们称为龙贝格积分法。
计算方法--龙贝格积分法
结论:从运行结果可见， 当设定误差 eps=1e-5 时， 为了满足误差要求， 运行 n=14 时才能满足， 这样精度比较高，但是代价就是运行的时候速度有点慢。
按照题目要求， 要求尽可能的使结果精确一些， 因此我这里采用的是根据误差要求进行的计 算，计算程序如下： function [R,T,n]=Romberg3(myfunc,a,b,eps) %This function uses Romberg to calculate integration %myfunc is the integrand,a is integral lower bound %b is integral upperbound %eps is tolerance format long; n=5; while 1 T(n,:)=0; T(1,1)=(b-a)*(myfunc(a)+myfunc(b))/2; for k=1:n-1 sum=0; for i=1:2^k sum=sum+myfunc(a+(2*i-1)*(b-a)/(2^(k+1))); end T(k+1,1)=0.5*T(k,1)+(b-a)*sum/(2^(k+1)); end for i=2:n T(i,2)=T(i,1)+1/3*(T(i,1)-T(i-1,1)); end for i=3:n T(i,3)=T(i,2)+1/15*(T(i,2)-T(i-1,2)); end for i=4:n T(i,4)=T(i,3)+1/(4^3-1)*(T(i,3)-T(i-1,3)); end if abs(T(n,4)-T(n-1,4))<eps break; else n=n+1; end end R=T(n,4); 以上便是这次编写的龙贝格积分的程序。 在matlab的command window运行如下： >> myfuc=inline('log(1+x)/(1+x^2)','x'); >> eps=1e-5;
讲：ch6-4(龙贝格求积)知识分享
龙贝格算法
用若干个积分近似值推算出更为精确的 积分近似值的方法，称为外推方法。
序 列 TN,SN,CN,RN分 别 称 为 梯 形 序 列
辛 普 森 序 列 、列柯和特龙斯贝序格 序 列 。
由龙贝格序列还可续以外继推，但因
4m 4m 1
1及4m11
0(m
4),故通常
0.9460835
T(1) 2
16T1(2) T1(1) 15
0.9460832
例题3
T(0) 3
64T2(1) T2(0) 63
0.9460832
T 3 (0)T 2 (0)
210 7110 6 2
1 sinx
0
x
dx0.9460832
算法
f[x_]:=Sin[x]/x; a=0;b=1;f[0]=1;f[1]=0.8414709; h[k_]:=(b-a)/2^k;T[0,0]=(b-a)/2*(f[a]+f[b]); T1][+0_h,k[k_]]*:=S1u/m2[*fT[[a0+,kh-[k]*(2i-1)],{i,1,2^(k-1)}]; N[Table[T[0,k],{k,0,3}],6]; MatrixForm[%] T1,[km])_/,k(4_^]:=m(-41^);m*T[m-1,k+1]-T[mN[Table[T[m,k],{m,1,3},{k,0,2}],7]; MatrixForm[%];N[Table[T[1,k],{k,0,2}],7]; MatrixForm[%];N[Table[T[2,k],{k,0,1}],7]; MatrixForm[%];N[T[3,0],7]
梯形法的逐次分半算法
《实验5：龙贝格积分》
实验名称实验4实验地点6A-XXX 实验类型设计实验学时 2 实验日期20 /X/X ★撰写注意：版面格式已设置好（不得更改），填入内容即可。
一、实验目的
1.龙贝格积分
二、实验内容
1.实验任务
1.龙贝格积分
2.程序设计
1）数据输入（输入哪些数据、个数、类型、来源、输入方式）
const double eps=0.000005;
const int MaxD=100;
double f(double), double a, double b
2）数据存储（输入数据在内存中的存储）
函数
double f1(double x)
double f2(double x)
double f3(double x)
void Romberg(double f(double), double a, double b)
3）数据处理（说明处理步骤。
若不是非常简单，需要绘制流程图）
1.输入要处理的数据进入变量中
2.进行函数处理
3.输出函数处理结果
4）数据输出（贴图：程序运行结果截图。
图幅大小适当，不能太大）。
计算方法-龙贝格资料
f
(2n )
若f ( x)在[a, b]变 化 不 大 ， 即f (n ) f (2n )， 则
b
a
b
f ( x)dx Tn
4
a f ( x)dx T2n
计算方法 © 2009, Henan Polytechnic University
5 July 2020
2
第六章 数值积分与数值微分
从而
b a
f
( x)dx
T2n
13（T2n
Tn）
即 ：b a
f
( x)dx
T2n
13（T2n
Tn）
由此可以认为，当T2n Tn 时,
b
a f ( x) T2n
计算方法 © 2009, Henan Polytechnic University
5 July 2020
3
第六章 数值积分与数值微分
5 July 2020
4
第六章 数值积分与数值微分
设h b a ,则 n
x0 x1 x2 x3
xn2 xn1 xn
T2n
h 4
n 1
{ f (a) 2 [
k 1
n 1
f (xk )
k 0
f
(
x k
1
)]
2
f (b)}
h 4
f
(a)
2
n1 k 1
f ( xk )
f (b)
h 2
计算方法 © 2009, Henan Polytechnic University
5 July 2020
1
第六章 数值积分与数值微分
若
用Tn及T2
分
《龙贝格求积》课件
龙贝格求积法的历史背景
01
龙贝格求积法是由瑞典数学 家龙贝格在19世纪末提出的
。
02
它的出现为数值积分的发展 奠定了基础，成为数学领域
的一项重要成果。
03
随着计算机技术的发展，龙 贝格求积法在科学计算、工 程技术和数据分析等领域得
缺点
计算量大
对于大规模问题，龙贝格求积法可 能需要较大的计算资源和时间，因 为需要进行多次迭代和数值逼近。
对初值敏感
该方法对初值的选择比较敏感，如 果初值选择不当，可能会导致算法
不收敛或者收敛到非期望的解。
需要选择合适的参数
龙贝格求积法的精度和稳定性与参 数的选择密切相关，需要仔细选择 合适的参数值。
通过差商的方式计算插值多项式的导数。
导数在数值分析中的应用
用于估计函数的局部变化趋势，提高数值计算的精度。
龙贝格求积公式的推导
03
龙贝格求积公式的定义
龙贝格求积公式的推导过程
龙贝格求积公式的应用
利用插值多项式及其导数构造求积公式， 用于数值积分。
通过构造插值多项式和其导数，利用牛顿 -莱布尼茨公式推导得出。
对未来研究的展望
01
入 ，未来可以对龙贝格求积法进 行进一步的改进和优化，以提 高其计算效率和适用范围。
02
与其他方法的比较研究
可以开展龙贝格求积法与其他 积分方法的比较研究，深入探 讨各种方法的优缺点和应用场 景，为实际应用提供更加全面
的理论支持。
03
扩展应用领域
对不规则区域处理困难
对于不规则积分区域，龙贝格求积 法可能需要进行额外的处理和调整 ，增加了计算的复杂度。
龙贝格求 积分
龙贝格（Romberg ）求积法1.算法理论Romberg 求积方法是以复化梯形公式为基础，应用Richardson 外推法导出的数值求积方法。
由复化梯形公式 )]()(2)([2222b f h a f a f h T +++=可以化为)]()]()([2[212112h a f h b f a f hT +++=＝)]([21211h a f h T ++一般地，把区间[a,b]逐次分半k －1次，（k ＝1,2,……,n ）区间长度（步长）为kk m a b h -=，其中mk ＝2k －1。
记k T ＝)1(k T 由)1(k T ＝]))12(([21211)1(1∑=---++km j k k k h j a f h T 从而⎰badxx f )(＝)1(kT－)(''122k f h a b ξ- （1）按Richardson 外推思想，可将（1）看成关于k h ，误差为)(2k h O 的一个近似公式，因而，复化梯形公式的误差公式为⎰badxx f )(－)1(k T ＝......4221++k k h K h K ＝∑∞=12i i k i h K （2）取1+k h ＝k h 21有 ⎰ba dx x f )(－)1(1+k T ＝∑∞=+121221i ik ii hK （3）误差为)(2jh O 的误差公式 )(j kT＝)1(-j kT＋141)1(1)1(------j j k j k T T2.误差及收敛性分析（1）误差，对复化梯形公式误差估计时，是估计出每个子区间上的误差，然后将n 个子区间上的误差相加作为整个积分区间上的误差。
（2）收敛性，记h x i =∆，由于∑=++=ni i i n x f x f h f T 01))]()([2)(＝))()((21101∑∑-==∆+∆n i ni i i i i x x f x x f上面两个累加式都是积分和，由于)(x f 在区间],[b a 上可积可知，只要],[b a 的分划的最大子区间的长度0→λ时，也即∞→n 时，它们的极限都等于积分值)(f I 。
数值计算方法 龙贝格求积公式 - 龙贝格求积公式
类似地可验证：
3
S2
4
4
1
T4
4
1
1
T2
龙
贝 格 算 法
Sn
4
4
1
T2
n
4
1
1
Tn
即
Sn
4 3 T2n
1 3
Tn
Sn __ 辛 普 森 积 分 值
注 意 复 化 辛 普 森 求 积 公式 的 余 项
3
Rn
I
Sn
ba 180
h 2
4
f
(4) ( )
O(h4
)
龙 贝 格
R2 n
I
S2n
T4
1 2
T2
1 4
[
f
(
1 4
)
f
(
3 4
)]
0.9445135
典型例题
例2
用复化梯形公式计算积分I 1 dx
0 (1 x ) x
精确至三位有效数字
I
1 dx 0 (1 x )
x t x
1 2dt 0 1 t2
1
g(t )dt
0
T1
1 [g(0) 2
g(1)] 1.5
T2
1 2
[T1
1
g(0.5)] 1.55
T4
1 2
[T2
1 ( g(0.25) 2
g(0.75))] 1.5656
典型例题
T8
1 2 [T4
1 ( g(0.125) 4
g(0.375)
g(0.625)
g(0.875))]
1.5695
注意到:
T8 T4
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
在等距基点的情况下，用计算机计算积分值通常都采用把区间逐次分半的方法进行。这样，前一次分割得到的函数值在分半以后仍可被利用,且易于编程。
一．程序框图
二．源码
1.
f=inline('1/(x+1)');%输入函数
a=0;b=1;%取值边界
eps=10^(-7);
h=b-a;
R(1,1)=h*(f(a)+f(b))/2;
j=0;
err=1;
m=1;
while err>eps
j=j+1;
h=h/2;
S=0;
for i=1:m
x=a+h*(2*i-1);
a=0;b=pi/2;%取值边界
eps=10^(-7);
h=b-a;
R(1,1)=h*(2/pi)/2;
j=0;
err=1;
m=1;
while err>eps
j=j+1;
h=h/2;
S=0;
for i=1:m
x=a+h*(2*i-1);
S=S+f(x);
end
m=2*m;
R(j+1,1)=R(j,1)/2+h*S;
龙贝格方法是在积分区间逐次二分的过程中，通过对梯形之值进行加速处理，从而获得高精度的积分值。
准备初值 和 ，用梯形计算公式计算出积分近似值
按区间逐次分半计算梯形公式的积分近似值,令
，
计算
，
按下面的公式积分（为便于编程，写下列形式）
梯形公式：
辛普生公式：
龙贝格公式：
步骤4 精度控制
当 ，（ 为精度）时，终止计算，并取 为近似值
3.
f=inline('log(x+1)/x');%输入函数
a=0;b=1;%取值边界
eps=10^(-7);
h=b-a;
R(1,1)=h*(1+log(2))/2;
j=0;
err=1;
m=1;
while err>eps
j=j+1;
h=h/2;
S=0;
for i=1:m
x=a+h*(2*i-1);
S=S+f(x);
S=S+f(x);
end
m=2*m;
R(j+1,1)=R(j,1)/2+h*S;
for i=1:j
R(j+1,i+1)=R(j+1,i)+(R(j+1,i)-R(j,i))/(4^i-1);
end
err=abs(R(j+1,j)-R(j+1,j+1));
end
ans=vpa(R(j+1,j+1),7)
for i=1:j
R(j+1,i+1)=R(j+1,i)+(R(j+1,i)-R(j,i))/(4^i-1);
end
err=abs(R(j+1,j)-R(j+1,j+1));
end
ans=vpa(R(j+1,j+1),8)
三．结果
S=S+f(x);
end
m=2*m;
R(j+1,1)=R(j,1)/2+h*S;
for i=1:j
R(j+1,i+1)=R(j+1,i)+(R(j+1,i)-R(j,i))/(4^i-1);
end
err=abs(R(j+1,j)-R(j+1,j+1));
end
ans=vpa(R(j+1,j+1),7)
2.
f=inline('log(x+1)/(x^2+1)');%输入函数
a=0;b=1;%取值边界
eps=10^(-7);
h=b-a;
R(1,1)=h*(f(a)+f(b))/2;
j=0;
err=1;
m=1;
while err>eps
j=j+1;
h=h/2;
S=0;
for i=1:m
x=a+h*(2*i-1);
end
m=2*m;
R(j+1,1)=R(j,1)/2+h*S;
for i=1:j
R(j+1,i+1)=R(j+1,i)+(R(j+1,i)-R(j,i))/(4^i-1);
end
err=abs(R(j+1,j)-R(j+1,j+1));
end
ans=vpa(R(j+1,j+1),7)
4.
f=inline('sin(x)/x');%输入函数