Romberg算法
龙贝格算法
龙贝格积分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 =+-=-应当是更好的结果。
龙贝格算法
2）理查森外推加速 从梯形公式出发， 将区间[a,b]逐次二分可以提高求积公式精度， 当[a,b] 分为 n 等份时，若记T������ = T ℎ ，当区间[a,b]分为 2n 等份时，则有 T2������ = T
ℎ 2
。再有泰勒公式展开为： T ℎ = ������ + ������1 ℎ2 + ������2 ℎ+ ⋯ +
建立一个命名为Romberg.m的function文件：
function[T]=Romberg(f,a,b,e) T=zeros(10,10); T(1,1)=(b-a)/2*(f(a)+f(b)); for k=2:10, sum=0; for i=1:2^(k-2), x=a+(2*i-1)*(b-a)/2^(k-1); sum=sum+f(x); end T(k,1)=T(k-1,1)/2+[(b-a)/2^(k-1)]*sum; for j=2:k, end %第一列用递推梯形公式 %定义龙贝格函数 %定义10阶的零元矩阵
1 3/2 ������ 0
������������。
算，并取������ ������, ������ ≈ ������ ；否则令 k=k+1，转（3）继续计算。 6）下图为我按照自己的算法所设计的示意表： 算法设计表：
k 1 2 3 4 …
h b-a (b-a)/2 (b-a)/4 (b-a)/8 …
������ ������, 1 T(1,1) T(2,1) T(3,1) T(4,1)
������ ������, 2
������ ������, 3
������ ������, 4
������ ������, 5
5.3龙贝格(Romberg)算法
一、复合梯形公式的递推化
将定积分I = ∫ f ( x )dx的积分区间[ a , b ]分割为n等份 a b−a xk = a + jh , j = 0 ,1,L , n 各节点为 h= n 复合梯形(Trapz)公式为
n−1 b−a Tn = [ f ( a ) + 2 ∑ f ( x j ) + f (b )] 2n j =1
k −3
if(fabs(t3[k-3]-t2[k-2])<0.000000001) // 若 | R2 { printf("k=%3d,I=%12.9f\n",k,t3[k-3]); break; } } } double f(double x) { double y; y=sqrt(1.0+x*x); return y; }
13
其中外推加速公式可简化为
1 Tm ( k − 1) = m [ 4 m Tm − 1 ( k ) − Tm − 1 ( k − 1)] 4 −1
--------(11)
并且m可以推广到 m = 1,2 ,L
Romberg算法求解步骤
k = 1 , 2 ,L
Romberg算法的代 数精度为m的两倍 Romberg算法的收敛 阶高达m+1的两倍
第6章 数值积分与数值微分
龙贝格(Romberg)算法 算法 龙贝格
1
龙贝格(Romberg)算法 算法 龙贝格
综合前几节的内容,我们知道 梯形公式,Simpson公式,Cotes公式的代数精度分别为 1次,3次和5次 复合梯形、复合Simpson、复合Cotes公式的收敛阶分别为 2阶、4阶和6阶 无论从代数精度还是收敛速度,复合梯形公式都是较差的 有没有办法改善梯形公式呢?
7.3Romberg积分(共64张PPT)
第五页，共六十四页。
表7-1
O(h)
O(h2 )
O(h3 )
O(h4 )
例1 设
带余项的差分(chà fēn)公式为
第六页，共六十四页。
导出具有误差为 解令
用 h/2代替(dàitì)h,得
(11) 的外推公式.
h2
(12)
为消去含 h2的项,用4乘(12)式减去 (11)式,得
第七页，共六十四页。
第二页，共六十四页。
把(1)式改写(gǎixiě)为 (3)
用h/2代替(3)式中的h,得
(4) 用2乘(4)式再减去(3)式,消去含h的项,得
(5)
令
,且记
第三页，共六十四页。
那么(5)式可写为
(6) 这里, 逼近 的误差为
再用 h/2 代替(dàitì) h , 使(6)式变为
(7) 用4乘(7)式减去(6)式,消去含 的项,得
则所谓 gn (x)与 的根相同,即是指这两个正
交多项式的根有如下的关系.
第三十二页，共六十四页。
性质5 (1) 区间[a,b]上带权函数 (的x) 正交
多项式序列
与
对应元素之间
只相差一个比例常数.
(2)区间[a,b]上带权函数
首项 系数 (x)
(shǒu xiànɡ)
为1的
正交多项式序列
唯一.
常见的正交多项式有Legendre(勒让德)多
7.5.1 引言 求积公式
当求积系数
、求积节点
(1) 都可以
自由选取时,其代数精确度最高可以达到多少 次?
下面的引理可以回答上述(shàngshù)问题.
第二十四页，共六十四页。
引理1 当求积系数
Romberg龙贝格算法实验报告
Romberg龙贝格算法实验报告课程实验报告课程名称：专业班级： CS1306班学号： Uxx14967 姓名：段沛云指导教师：报告日期：计算机科学与技术学院目录1 实验目的 ........................................................12 实验原理 ........................................................13 算法设计与流程框图 (2)4 源程序 ........................................................ .. 45 程序运行 ........................................................76 结果分析 ........................................................77 实验体会 ........................................................71 实验目的掌握Romberg公式的用法，适用范围及精度，熟悉Romberg算法的流程，并能够设计算法计算积分31得到结果并输出。
1x2 实验原理2.1 取k=0,h=b-a,求T0=数)。
2.2 求梯形值T0(b-a),即按递推公式（4.1）计算T0。
k2h[f(a)+f(b)]，令1→k,(k记区间[a，b]的二分次22.3 求加速值，按公式(4.12)逐个求出T表的第k行其余各元素Tj(k-j)(j=1,2,….k)。
2.4 若|Tk+1-Tk|n-111T2n=[Tn+hn∑f(xi+)]22i=01Sn=T2n+(T2n-Tn)31Cn=S2n+(S2n-Sn)151Rn=C2n+(C2n-Cn)633 算法设计与流程框图算法设计：(先假定所求积分二分最大次数次数为20) 3.1 先求T[k][0] 3.2 再由公式T(k)m4m(k+1)1)=mTm-1-mTm(k-1(k=1,2,) 求T[i][j] 4-14-13.3 在求出的同时比较T[k][k]与T[k-1][k-1]的大小，如果二者之差的绝对值小于1e-5，就停止求T[k][k]；此时的k就是所求的二分次数，而此时的T[k][k]就是最终的结果 3.4 打印出所有的T[i][j]；程序流程图4 源程序#include #include #include #include int main(void) {float f(float(x)) {float y; y=1/x; return y; }floata,b,e,h,s,k,x,T1=0,T2=0,S1=0,S2=0,C1=0,C2=0,R1=0,R2=0; inti=0;printf("请输入积分下限 : "); scanf("%f",&a);printf("\n请输入积分上限 :"); scanf("%f",&b);printf("\n请输入允许误差 :"); scanf("%f",&e); k大学网=1; h=b-a;T1=h*(f(a)+f(b))/2;printf("____________________________________________\n"); printf("计算结果如下 : \n");printf("\nk T2 S2 C2 R2\n");printf("%d %10.7f %10.7f %10.7f %10.7f\n",i,T1,S1,C1,R1); do {x=a+h/2; s=0; while(x{ s=s+f(x); x=x+h; }T2=(T1+s*h)/2; S2=T2+(T2-T1)/3; if(k==1) {T1=T2; S1=S2; h=h/2; k=k+1; }else if(k==2) {C2=S2+(S2-S1)/15; C1=C2; T1=T2; S1=S2; h=h/2; k=k+1; }else if(k==3) {R2=C2+(C2-C1)/63; C2=S2+(S2-S1)/15; C1=C2; T1=T2; S1=S2; h=h/2; k=k+1; } else {C2=S2+(S2-S1)/15;R2=C2+(C2-C1)/63; if(fabs(R2-R1)printf("%d %10.7f %10.7f %10.7f %10.7f\n",i+1,T2,S2,C2,R2); break;} else { R1=R2; C1=C2; T1=T2; S1=S2; h=h/2; k=k+1; } } i++; printf("%d %10.7f %10.7f %10.7f %10.7f\n",i,T2,S2,C2,R2); } while(1); system("pause"); return 0; }5 程序运行6 结果分析如上所示的结果与课本中求得的结果完全一样，表明程序编写正确，且符合要求，事实上，只要再将所求值的精度设置得更小，则所求的结果将更加准确，最终将无限接近于标准值，由上表也可以看出用龙贝格积分法求函数的积分值在精度比较低的情况下就能求到很准确的值！7 实验体会本次实验较为简单，主要时间是耗费在循环判断上面，因为书上已经给了流程图，都是基本的C语言，难度不大。
数学实验题目2 Romberg积分法
数学实验题目2 Romberg 积分法摘要考虑积分()()b aI f f x dx =⎰欲求其近似值，可以采用如下公式：（复化）梯形公式 110[()()]2n i i 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 →∞=这是一个收敛速度更快的一个数值求积公式，我们称为龙贝格积分法。
最小二乘法和Romberg求积分
计算方法实验报告湖北师范学院计算机科学与技术学院一、实验题目最小二乘法和Romberg求积分二、实验内容1、最小二乘法拟合本实验是对于已知的m+1对不带权系数{}mi iw=的离散数据{},mi i ix y=，计00,min()max()i m i ma bi ix y≤≤≤≤==。
在连续的函数空间[]()(),0C a b f a f b⨯<中选定n+1个线性无关的基函数(){}nk kxσ=来构造拟合函数()()nk k kkx a xσσ==∑，求得当系数是最优化模型()()()201min I,,,,mn k i iia a a x y xσσ=⋅⋅⋅=-⎡⎤⎣⎦∑时的解***01,,,na a a⋅⋅⋅。
2、romberg求积分法Romberg算法是利用复合梯形公式，在对积分区间的步长逐次折半的过程中，求得积分I=∫a b f(x)dx的近似值序列{T2k}。
三、法思想描述1、最小二乘法拟合（1）、根据最小二乘法的拟合原理将问题简化为求由待拟合离散数据的正则方程组的求解问题；（2）、输入一散数据{},mi i ix y=，存放于数组[],[]X m Y m中，并给定拟合次数n根据待拟合的次数确定对应的系数矩阵为范德蒙行列式的线性方程组；（3）、由[][][][]TA m n A m n⨯得到正则方程组的系数矩阵[][]F n n，如下式，同时由[][][]TA m n Y m⨯得到值向量[]B n；(4)、用最列主元素法解出正则方程组的根，就得到“二”中描述的系数的最优解。
2、romberg求积分法Romberg算法是区间主次二分的过程中，对梯形值进行加权平均以获得准确程度较高a0a1a2…a n1x0x0 2…x0n1 x1x i n…x i n……………………1 x m x m2…x m n＝Y1Y2LLLY m的积分值的一种方法。
对于定积分I= 的复化梯形公式，其余项将积分区间[a，b]，逐次折半，假设,以保证复化梯形公式余项系数是非零的，则构成相应的外推法称为Romberg算法：下标m外推得到的第m个算法。
龙贝格(Romberg)算法的应用实验报告
Lab4 龙贝格（Romberg）算法的应用下面图1中的塑料雨蓬材料是由图2中所示的长方形平板塑料材料压制而成。
图1 图2已知图1的横截面曲线形状满足函数，则给定了雨蓬的长度后，要求需要平板原材料的长度。
函数接口定义：double Integral(double a, double b, double (*f)(double x, double y, double z), double TOL, double l, double t)在接口定义中：a、b分别为定积分的上、下界，f是积分核函数，其中x是积分哑元，y、z是本题目定义的特殊参数，分别对应中的l和t；TOL是要求积分达到的精度；l和t传入裁判输入的参数l和t的值。
另注意：的单位是厘米，输出的长度值要求以米为单位。
裁判程序样例如下：#include<stdio.h>#include<math.h>double f0( double x, double l, double t ){ /* 弧长积分的核函数*/return sqrt(1.0+l*l*t*t*cos(t*x)*cos(t*x));}double Integral(double a, double b, double (*f)(double x, double y, double z), double TOL, double l, double t);int main(){double a=0.0, b, TOL=0.005, l, t;while (scanf("%lf %lf %lf", &l, &b, &t) != EOF)printf("%.2f\n", Integral(a, b, f0, TOL, l, t));return 0;}裁判输入样例：2 100 1标准输出样例：1.68实验报告：1.求解步骤参照书上的龙贝格求积算法2.步骤①利用k=0;h=b-a;T[0][0]=(h/2)*(f(a,l,t)+f(b,l,t));求解T0(0)3.求解T0(0)到T0(k)的值由公式h=b−an 及h=b−a2k可得n=2k其中h为步长，n为二分次数又由递推公式：T2n=12T n+ℎ2∑f(xk+12)n−1k=0得T2k+1=12T2k+b−a2k+1∑f[a+b−a2k+1(2i−1)]2k−1i=1，k=0,1,2,3~其中xk+12= a+ℎ2(2i−1)。
Romberg算法
上例说明{Tn}收敛慢，求T128 要计算64个新增的函 数值，而将T8与T4重新组合可构造S8。
1 4 1 T23 (T23 T2 2 ) T8 T4 S 4 3 3 3
K
T2k
S2k-1 0.9461459 0.9460869 0.9460833 …
m
C2k-2
R2k-3
因此计算梯形序列{T2m}可按：
T2m hm ( f (a) f (b) 2 f (a khm )) 2 k 1
2m1
注意到每个子区间 xk , xk 1 经过等分只增加了一个分点： 1 x 1 ( xk xk 1 ) 用复化梯形公式可求得 xk , xk 1 上的积分值为 k 2 2
1 ~ 4 如, 当n 1 时， T T2 T1 3 3 4ba ab 1ba [ f (a) 2 f ( ) f ( b )] [ f ( a ) f ( b )] 3 4 2 3 2 ba ab [ f (a) 4 f ( ) f ( b )]. 6 2 4 1 即 S1 T2 T1 . 3 3
3、加密一次区间 ( 再分半)h2
ba ba , xk a kh2 a k, k 0,1, 2, 3,4 4 4 3 3 ba 1 1 h2 T4 ( f ( a ) ( f ( a kh2 ) f ( b )) ( f ( a ) f ( b ) 2 f ( a kh2 ) k 1 k 1 4 2 2 2
1 h n1 T2n Tn f ( xk 1 ) 2 2 2 k 0
此为复化梯形公式的递推公式
例2 利用变步长的梯形法求I dx 的近似值. x 1 T1 [ f (0) f (1)] 0.9207355 . 2 1 1 1 T2 T1 f ( ) 0.9397933 . 2 2 2 1 1 1 3 T4 T2 [ f ( ) f ( )] 0.9445135 . 2 4 4 4 1 1 1 3 5 7 T8 T4 [ f ( ) f ( ) f ( ) f ( )] 0.956909 . 2 8 8 8 8 8
利用Matlab实现Romberg数值积分算法----系统建模与仿真结课作业
利用Matlab 实现Romberg 数值积分算法一、内容摘要针对于某些多项式积分，利用Newton —Leibniz 积分公式求解时有困难，可以采用数值积分的方法，求解指定精度的近似解，本文利用Matlab 中的.m 文件编写了复化梯形公式与Romberg 的数值积分算法的程序，求解多项式的数值积分，比较两者的收敛速度。
二、数值积分公式1.复化梯形公式求解数值积分的基础是将区间一等分时的Newton —Cotes求积公式:I =(x)[f(a)f(b)]2bab af dx -≈+⎰ 其几何意义是，利用区间端点的函数值、与端点构成的梯形面积来近似(x)f 在区间[a,b]上的积分值，截断误差为：3"(b a)()12f η-- (a,b)η∈ 具有一次的代数精度，很明显，这样的近似求解精度很难满足计算的要求，因而，可以采用将积分区间不停地对分，当区间足够小的时候，利用梯形公式求解每一个小区间的积分近似值，然后将所有的区间加起来，作为被求函数的积分，可以根据计算精度的要求，划分对分的区间个数，得到复化梯形公式：I =11(b a)(b a)(x)dx [f(a)f(b)2(a )]2n bak k f f n n -=--≈+++∑⎰其截断误差为：2"(b a)h ()12R f η--=(a,b)η∈ 数值积分算法使用复化的梯形公式计算的数值积分，其收敛速度比减慢，为此，采用Romberg 数值积分。
其思想主要是，根据I 的近似值2n T 加上I 与2n T 的近似误差，作为新的I 的近视，反复迭代，求出满足计算精度的近似解。
用2n T 近似I 所产生的误差可用下式进行估算：12221()3n n n I T T T -∆=-=-新的I 的近似值：122n n j T T -=∆+ j =（0 1 2 ….） Romberg 数值积分算法计算顺序i=0 (1) 002Ti=1 (2) 102T(3) 012Ti=2 (4) 202T (5) 112T(6) 022Ti=3(7)302T (8) 212T (9) 122T(10) 032T i=4 (11)402T(12) 312T(13) 222T(14) 132T…………其中，第一列是二阶收敛的，第二列是四阶收敛的，第三列是六阶收敛的，第四列是八阶收敛的，即Romberg 序列。
Romberg算法
Romberg算法张冰 229201322039231．Romberg算法程序%f是被积函数，a是积分上限，b是积分下限，e是精度function Romberg(f,a,b,e)format longm=100;T=zeros(m,m); %初始化一个100*100的矩阵，当精度满足要求时会确定实际的矩阵大小T(1,1)=(b-a)*(f(a)+f(b))/2; %计算T(1,1)for i=1:msum=0;for j=1:2^(i-1)sum=sum+f(a+(2*j-1)*(b-a)/(2^i));endT(i+1,1)=0.5*T(i,1)+(b-a)*sum/2^i; %计算T(2,1)for k=1:iT(i+1,k+1)=(4^k*T(i+1,k)-T(i,k))/(4^k-1); %外推其他T(i,k) endif abs(T(i+1,i+1)-T(i+1,i))<=e %精度满足要求时停止计算break;endendfor i=1:mif T(i,1)==0 %当第一列第i个元素为0时，矩阵大小确定为ij=i;break;endenddisp('计算表格如下：')T=T(1:j-1,1:j-1) %输出计算结果矩阵disp('所求的积分值为：')disp(T(j-1,j-1)) %输出最终计算值PS：这个程序写起来还是比较容易的，基本上就是把书上给的龙贝格矩阵计算方法用MATLAB语言来实现就可以了，但是完成的时候遇到两个问题：第一个是因为矩阵下表非常多，所以容易出错导致计算结果有误；第二个是最开始不知道矩阵的实际大小是多少，我想到的办法就是先初始化一个足够大的矩阵，虽然觉得这可能不是一个好办法……2．利用Romberg完成书上的例6%输入要计算的函数function y = f(x)y=x^(3/2);end计算结果如下：Ps：关于这个输出结果我有几点觉得不太满意但是改了很久都不能实现……第一个是输出矩阵每列里的0不能对齐，所以看着会很别扭，但是又不知道怎么不输出0；第二个是书上给的例子中保留了六位小数，但是MATLAB我只会保留四位小数或者十五位小数，所以书上的最终结果是0.400002，我得出的结果是0.400001516355028……。
romberg算法
romberg算法
普朗伯格算法（Pratt-Romburg Algorithm）是一种用于解决最小二乘问题的数值计算方法。
它是由美国数学家Ralph E. Pratt和美国数学家John Romburg于1962年提出的。
普朗伯格算法的基本思想是，将最小二乘问题转化为一个最优化问题，然后使用梯度下降
法来求解。
梯度下降法是一种迭代算法，它通过不断更新参数来最小化损失函数。
普朗伯格算法的优点是，它可以解决非线性最小二乘问题，而且它的计算复杂度比其他算
法要低。
它的缺点是，它的收敛速度比其他算法要慢，而且它可能会收敛到局部最优解。
普朗伯格算法在机器学习领域有着广泛的应用，它可以用于拟合模型，估计参数，求解最优化问题等。
它也可以用于解决统计学中的最小二乘问题，如线性回归、非线性回归等。
总之，普朗伯格算法是一种有效的数值计算方法，它可以用于解决最小二乘问题，并且在机器学习和统计学领域有着广泛的应用。
数值分析中的Romberg积分误差控制技巧研究方向
数值分析中的Romberg积分误差控制技巧研究方向数值分析中的Romberg积分误差控制技巧研究方向在数值分析领域中，Romberg积分是一种重要的数值积分方法，常用于解决复杂函数的数值积分问题。
然而，由于计算机运算的误差以及数值积分方法本身的近似性质，Romberg积分在实际应用中的误差控制成为一个需要深入研究的问题。
本文将围绕Romberg积分误差控制技巧展开讨论。
1. Romberg积分方法简介Romberg积分方法是一种基于复化梯形法的数值积分方法，在求解定积分时具有较高的精度。
其基本思想是通过不断提高复化的精细程度，逼近定积分的真实值。
Romberg积分方法的核心是构造一个误差序列，通过这个序列的收敛性来控制数值积分的误差。
2. Romberg积分误差控制的方法2.1 自适应步长调整在Romberg积分中，通过将积分区间等分为一系列子区间，并分别计算每个子区间的积分值，然后通过合并这些子区间的积分值来逼近定积分。
自适应步长调整的方法是在每次迭代计算中，根据当前计算的误差估计值来自动调整子区间的精细程度。
通过这种方式，可以在保证计算效率的同时，控制积分误差的大小。
2.2 Richardson外推Richardson外推是一种常用的误差控制技巧，可以在Romberg积分中得到广泛应用。
该方法通过利用不同精度的近似值之间的误差差异来估计积分的误差，并通过外推计算来提高积分的精度。
Richardson外推方法的关键是选择合适的外推系数，使得误差得到更精确的估计。
2.3 自适应不等距节点选取不等距节点选取是Romberg积分中另一种常用的误差控制技巧。
通常情况下，积分节点以等距方式选取，但这种方式可能导致积分误差较大。
因此，自适应不等距节点选取方法通过优化节点的位置，使得在积分区间上更加均匀地分布积分点，从而提高积分的精度。
3. 数值实验与比较为了验证上述误差控制技巧的有效性，我们进行了一系列数值实验并对结果进行了比较。
44Romberg积分法
一、Richadson外推加速法 二、Romberg公式 三、梯形法的递推化
一、李查逊(Richardson)外推加速法
复化梯形法的误差公式当h → 0时为：
I
Tn
ba 12
h2
f
( )
即积分值Tn的截断误差大致为O(h2)，事实上 有
定理 设 f ( x) C [a,b] ，T (h) Tn 则成立
，在每个小区间[xk，xk+1]上
用梯形公式计算为
Ik
h 2
fቤተ መጻሕፍቲ ባይዱ
(
xk
)
f ( xk1 )
如果再二分一次，则步长减半，即h/2，分点
增至2n+1个，记区间[xk，xk+1]上
经过二分后新增分点为
xk
1 2
1 2
(
xk
xk 1)
用复化梯形公式求得该区间上的积分值为
h/ 2
2
f
(
xk
)
2
f
(
xk
1 2
C
2n
1 63
C
n
综合上面的加工过程，有
对梯形公式加工： 对辛甫生公式加工： 对柯特斯公式加工：
Sn
4 3 T2n
1 3 Tn
16
1
Cn 15 S2n 15 Sn
Rn
64 63 C2n
1 63 Cn
实质：将粗糙的梯形公式值逐步加工成精度较 高的公式。
结论：
来源于：李查逊（Richardson） 外推加速法
定义：利用李查逊（Richardson）外推加速法将低 阶求积公式通过加密节点，再加工成高阶的求积公 式，统称为龙贝格（Romberg）公式。
Romberg求积分公式
《MATLAB程序设计实践》课程考核1、编程实现以下科学计算算法，并举一例应用之。
“Romberg求积分公式”2、编程解决以下科学计算和工程实际问题。
1）、给定半径的为r，重量为Q的均质圆柱，轴心的初始速度为v0，初始角速度为w0且v0>r*w0，地面的摩擦系数为f，问经过多少时间后，圆柱将无滑动地滚动，求此时的圆柱轴心的速度。
2）、在一丘陵地带测量高程，x和y方向每隔100m测一个点，得高程数据如下，试拟合一曲面确定合适的模型，并由此找出最高点和该点的高程。
一、Romberg求积分公式1、算法说明：此算法可自动改变积分步长,使其相临两个值的绝对误差或相对误差小于预先设定的允许误差.Romberg加速法公式在等距节点的情况下，通过对求积区间（a，b）的逐次分半，由梯形公式出可逐次提高求积公式精度，这就是Romberg求积的基本思路，由于梯形公式余项只有精度，即，但当节点加密时可组合成其精度达到，如果再由与组合成则可使误差精度达到，于是依赖于x，若在上各阶导数存在，将展开，可将展成的幂级数形式，即，记的计算精度，可利用外推原理逐次消去式右端只要将步长h逐次分半，利用及组合消去，重复同一过程最后可得到递推公式，此时.说明用其误差阶为，这里表示m次加速。
计算时用序列表示区间分半次数，即具体计算公式为，就是Romberg求积方法。
2、程序代码：M文件1）、Romberg加速法function [s,n]=rbg1(a,b,eps)if nargin<3,eps=1e-6;ends=10;s0=0;k=2;t(1,1)=(b-a)*(f(a)+f(b))/2;while (abs(s-s0)>eps)h=(b-a)/2^(k-1);w=0;if (h~=0)for i=1:(2^(k-1)-1)w=w+f(a+i*h);endt(k,1)=h*(f(a)/2+w+f(b)/2);for l=2: kfor i=1:(k-l+1)t(i,l)=(4^(l-1)*t(i+1,l-1)-t(i,l-1))/(4^(l-1)-1);endends=t(1,k);s0=(t(1,k-1));k=k+1;n=k;else s=s0;n=-k;endend2）、改进的Romberg求积函数function [s,eer]=rbg2(a,b,eps)if nargin<3,eps=1e-6;endm=1;t(1,1)=(b-a)*(f(a)+f(b))/2;r(1,1)=0;while ((abs(t(1,m)-r(1,m))/2)>eps)c=0;m=m+1;for j=1:2^(m-1)c=c+f(a+(j-0.5)*(b-a)/2^(m-1));endr(m,1)=(b-a)*c/2^(m-1);for j=2:mfor k=1:(m-j+1)r(k,j)=r(k+1,j-1)+(r(k+1,j-1)-r(k,j-1))/(4^(j-1)-1);endt(1,j)=r(1,j-1)+2*(4^(j-2)-1)*(t(1,j-1)-r(1,j-1))/(4^(j-1)-1);endenderr=abs(t(1,m)-r(1,m))/2;s=t(1,m);3）定义f.m函数如下：function f=f(x);f=x.^3;4）运行命令及结果>> rbg1(0,2)>> rbg2(0,2)3、流程图二、圆柱体问题1、问题分析圆柱体水平方向受到地面的摩擦阻力f*Q，该摩擦力对轴心速度起减速作用，同时又产生一个力矩，对角速度起加速作用。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
Richardson 外推算法
Romberg 算法
记: T
(k) 0
T2k , T
(k) 1
S2k , T
(k) 2
C2k , T
(k) 3
R2k
T0( k ) : k 次等分后梯形公式计算所得的近似值 ( Tmk ) : m 次加速后所得的近似值
(0)
(1) (2) (3)
( Tmk ) ( ( 4 mTmk 1) Tmk )1 1 4m 1
ba h n
h0 b a
h1 ba 2
n=1
n=2 n=4
1 h1 1 T4 T2 f (a ih1 0.5h1 ) 2 2 i 0 1 h2 3 T8 T4 f ( a ih2 0.5h2 ) 2 2 i 0
ba h2 4
梯形法递推公式
Cha4.2.m
程序演示
Cha41.m: 复化N-C公式 Cha42.m: 梯形法的递推
龙贝格 (Romberg) 加速
梯形法递推公式算法简单，编程方便
但收敛速度较 慢 解决方法：龙贝格 (Romberg) 加速
复化积分公式的渐近状态
思想：利用余项公式与积分的定义
梯形公式：
I [ f ]=0.946083070367
1
0
sin( x ) dx x
T4 0.944513522
k
0 1 2 3 4 5 6 7 8 9 10
T2 0.939793285
T0(k)
0.920735492 0.939793285 0.944513522 0.945690864 0.945985030 0.946058561 0.946076943 0.946081539 0.946082687 0.946082975 0.946083046
2h 6 I Cn [ f (5) (b) f (5) (a)] 945 46
Romberg加速
梯形公式的加速
I T2 n 1 I Tn 4
辛普森公式
辛普森公式的加速
I S2n 1 I S n 16 I C2 n 1 I Cn 64
I
4 1 4 1 而 S n T2 n Tn I T2 n Tn 3 3 3 3 T2 n Tn 事后估计法 I T2 n 3
数值分析及计算软件
Chap 4 数值积分与数值微分
4.1 Romberg 算法
Romberg 算法
利用复化梯形公式、复化Simpson公式、复化 Cotes公式等计算定积分时，如何选取步长 h 太 大: 太 小: 计算精度难以保证
？
计算量增大
解决办法：变步长算法 通常采取将区间不断对分的方法，即取 n = 2k ， 反复采用复化求积公式，直到所得到的计算结果 满足指定的精度为止。
① T1 =T0
② T2 =T0
③ S1 =T1
(0) (1) (2)
Romberg 算法是收敛的
⑥ C1 =T2
⑨ C2 =T2
(0) (1) (0)
④ T4 =T0
⑦ T8 =T0
⑤ S2 =T1
⑧ S4 =T1
⑩ R1 =T3
举例
例：计算定积分
T1 0.920735492
梯形法的余项展开式的推导
T (h) I[ f ] 1h2 2h4 i h2i
必要性: Richardson外推法的基础
推导方法：Taylor展开方法
作业
1.证明：梯形公式的Romberg加速为辛普森公式 2. 教材 P 104/11
hk 1 ba 2 k 1
举例
计算精度满足 | T2n Tn | 107
I [ f ]=0.946083070367
例：用梯形法的递推公式计算定积分 解：
T (0)
1
0
sin( x ) dx , 要求 x
ba f ( a) f ( b) 0.920735492 2 1 (0) h0 0 T (1) T f (a ih0 0.5h0 ) 0.939793285 2 2 i 0 1 (1) h1 1 ( 2) T T f (a ih1 0.5h1 ) 0.944513522 2 2 i 0 1 ( 2) h2 1 ( 3) T T f ( a ih2 0.5h2 ) 0.945690864 2 2 i 0
定义
如果一种复化求积公式 I n 满足下列关系
I In lim p C , C 0 h 0 h
则称该求积公式是 p 阶收敛的. 有如下误差估计式
h2 I Tn [ f ' (b) f ' (a)] 12 h4 I Sn [ f ' ' ' (b) f ' ' ' (a)] 4 180 2
T (h) I[ f ] 1h2 2h4 i h2i
ba h n
其中系数
k , k 1,2,3,...
与h无关
梯形法的加速
T (h) I[ f ] 1h2 2h4 i h2i I[ f ] O( h2 )
梯形法递推公式
将 [a, b] 分成 n 等分 [xi , xi+1] ，xi a i h,
h ba n n1 n1 h f ( x ) f ( x ) h f (a) 2 f ( x ) f (b) Tn i i i 1 2 2 i 0 i 1
步长折半：[xi , xi+1/2] ， [xi +1/2 , xi+1]
h T2 n f ( xi ) f ( xi 1 2 ) f ( xi 1 2 ) f ( xi 1 ) i 0 4 n1 h f ( xi ) 2 f ( xi 1 2 ) f ( xi 1 ) i 0 4 h n1 h n1 1 h n1 f ( xi ) f ( xi 1 ) f ( xi 1 2 ) Tn f ( xi 1 2 ) 4 i 0 2 i 0 2 2 i 0
1 S1 4T2 T1 3 0.946145882
1 S2 4T4 T2 3 0.946086934
1 C1 16 S2 S1 0.946083004 15
Cha43.m
举例
例：用 Romberg 算法计算定积分
满足
( ( ) | Tmk ) Tmk 1 | 107
I Tn 1 [ f ' (b) f ' (a)] 2 h 12
辛普森公式：
I Sn 1 [ f ' ' ' (b) f ' ' ' (a )] 4 4 h 180 2
Cotes公式：
I Cn 2 [ f (5) (b) f ( 5) (a)] h6 945 46
h h h h T I [ f ] 1 2 i 2 2 2 2
2 4 2i
4T ( h/2) T (h) 3I[ f ] ( 3 / 4)2h4 ( 15 ห้องสมุดไป่ตู้ 16)3h6
I [ f ]=0.4
1
0
x 3 dx , 要求计算精度
Cha44.m
解：逐步计算可得
k
0 1 2
T0(k)
0.50000000
T1(k)
T2(k)
T3(k)
T4(k)
T5(k)
0.42677670 0.40236893 0.40701811 0.40043192 0.40030278
3 3 3
0.40181246 0.40007725 0.40005361 0.40004965 0.40046340 0.40001371 0.40000948 0.40000878 0.40000862 0.40011767 0.40000243 0.40000168 0.40000155 0.40000152 0.40000152
hk 1 1 T2k T2k 1 2 2
2k 1 1
i 0
f ( a ihk 1 0.5hk 1 )
hk 1
ba k 1 2
k 1
记 T
(k)
T2k
2k 1 1
n2
T
(k)
1 ( k 1) hk 1 T 2 2
i 0
f (a ihk 1 0.5hk 1 )
n1
梯形法递推公式
1 h n1 1 h n1 T2 n Tn f ( xi 1 2 ) Tn f ( a ih 0.5h) 2 2 i 0 2 2 i 0
ba T1 f (a) f (b) 2 h0 0 1 T2 T1 f (a ih0 0.5h0 ) 2 2 i 0
Cotes公式