[计算机数值积分]龙贝格公式求数值积分

Spring-_-Bear 的 CSDN 博客导航

梯形法的算法简单，但精度低，收敛速度缓慢。如何提高收敛速度以节省计算量，自然是人们极为关心的课题。

根据梯形法的误差公式

I − T n ≈ − h 2 12 [ f ′ ( b ) − f ′ ( a ) ] I-T_{n}\approx -\frac{h^{2}}{12}[f'(b)-f'(a)] I−Tn≈−12h2[f′(b)−f′(a)]

积分值 T n T_{n} Tn 的截断误差大致与 h² 成正比，因此步长减半后误差将减至 1 4 \frac{1}{4} 41，即有

I − T 2 n I − T n ≈ 1 4 \frac{I-T_{2n}}{I-T_{n}}\approx \frac{1}{4} I−TnI−T2n≈41

将上式移项整理，知

I − T 2 n ≈ 1 3 ( T 2 n − T n ) (1) I-T_{2n}\approx \frac{1}{3}(T_{2n}-T_{n})\tag1 I−T2n≈31(T2n−Tn)(1)

由此可见，只要二分前后两个积分值 T n T_{n} Tn 与 T 2 n T_{2n} T2n 相当接近，就可以保证计算结果 T 2 n T_{2n} T2n 的误差很小。

按式（1），积分值 T 2 n T_{2n} T2n 的误差大致等于 1 3 ( T 2 n − T n ) \frac{1}{3}(T_{2n}-T{n}) 31(T2n−Tn)，如果用这个误差值作为 T 2 n T_{2n} T2n 的一种补偿，可以期望所得的

T ‾ = T 2 n + 1 3 ( T 2 n − T n ) = 4 3 T 2 n − 1 3 T n (2) \overline T=T_{2n}+\frac{1}{3}(T_{2n}-T_{n})=\frac{4}{3}T_{2n}-\frac{1}{3}T_{n}\tag2 T=T2n+31(T2n−Tn)=34T2n−31Tn(2)

应当是更好的结果。

再考察此例，所求得的两个梯形值 T 4 = 0.9445135 T_{4}=0.944 513 5 T4=0.9445135 和 T 8 = 0.9456909 T_{8}=0.945 690 9 T8=0.9456909 的精度都很差（与准确值 0.9460831 比较，它们只有二三位有效数字），但如果将它们按式（2）做线性组合，则新的近似值

T ‾ = 4 3 T 8 − 1 3 T 4 = 0.9460834 \overline T=\frac{4}{3}T_{8}-\frac{1}{3}T_{4}=0.946 083 4 T=34T8−31T4=0.9460834

却有六位有效数字。

按式（2）组合得到的近似值 T ‾ \overline T T，其实质究竟是什么呢？直接验证易知

S n = 4 3 T 2 n − 1 3 T n (3) S_{n}=\frac{4}{3}T_{2n}-\frac{1}{3}T_{n}\tag3 Sn=34T2n−31Tn(3)

这就是说，用梯形法二分前后两个积分值 T n T_{n} Tn 与 T 2 n T_{2n} T2n 按式（2）作线性组合，结果得到辛普生的积分值 S n S_{n} Sn。

再考察辛普生法。按式

I − S n ≈ − 1 180 ( h 2 ) 4 [ f ′ ′ ′ ( b ) − f ′ ′ ′ ( a ) ] I-S_{n}\approx -\frac{1}{180}(\frac{h}{2})^{4}[f'''(b)-f'''(a)] I−Sn≈−1801(2h)4[f′′′(b)−f′′′(a)]

其截断误差与 h 4 h_{4} h4 成正比。因此，若将步长折半，则误差相应地减少至 1/16，即有

I − S 2 n I − S n ≈ 1 16 \frac{I-S_{2n}}{I-S_{n}}\approx \frac{1}{16} I−SnI−S2n≈161

由此得

I ≈ 16 15 S 2 n − 1 15 S n I\approx \frac{16}{15}S_{2n}-\frac{1}{15}S_{n} I≈1516S2n−151Sn

不难验证，上式右端的值其实等于 C n C_{n} Cn，就是说，用辛普生法二分前后的两个积分值 S n S_{n} Sn 与 S 2 n S_{2n} S2n，按上式再作线性组合，结果得到柯特斯法的积分值 C n C_{n} Cn，即有

C n ≈ 16 15 S 2 n − 1 15 S n (4) C_{n}\approx \frac{16}{15}S_{2n}-\frac{1}{15}S_{n}\tag4 Cn≈1516S2n−151Sn(4)

重复同样的手续，依据柯特斯的误差公式

I − C n ≈ − 2 945 ( h 4 ) 6 [ f ( 5 ) ( b ) − f ( 5 ) ( a ) ] I-C_{n}\approx -\frac{2}{945}(\frac{h}{4})^{6}[f^{(5)}(b)-f^{(5)}(a)] I−Cn≈−9452(4h)6[f(5)(b)−f(5)(a)]

可进一步导出下列龙贝格公式

R n = 64 63 C 2 n − 1 63 C n (5) R_{n}=\frac{64}{63}C_{2n}-\frac{1}{63}C_{n}\tag5 Rn=6364C2n−631Cn(5)

应当注意，龙贝格公式（5）已经不属于牛顿 - 柯特斯公式的范畴。

在步长二分的过程中运行公式（3）~（5）加工三次，就能将粗糙的积分值 T n T_{n} Tn 逐步加工成精度更高的龙贝格值 R n R_{n} Rn，或者说，将收敛缓慢的梯形值序列 T n T_{n} Tn 加工成收敛迅速的龙贝格值序列 R n R_{n} Rn，这种加速算法称为龙贝格算法，其加工流程如图 2-3 所示：

龙贝格算法的程序框图如图 2-4 所示：

例：基于复化梯形的递推公式的变步长算法求积分中的例子，用龙贝格算法计算

I = ∫ 0 1 sin ⁡ x x d x I=\int_{0}^{1}\frac{\sin x}{x}dx I=∫01xsinxdx

得到的梯形值，计算结果如下表所示，表中 k 代表二分次数。

k	T₂^k	S₂^k-1	C₂^k-2	R₂^k-3
0	0.920 135 5
1	0.939 793 3	0.946 145 9
2	0.944 513 5	0.946 086 9	0.946 083 0
3	0.945 690 9	0.946 083 4	0.946 083 1	0.946 083 1

我们看到，这里用二分 3 次的数据（它们的精度都很低，只有二三位有效数字），通过三次加速获得了例中需要二分 10 次才能得到的结果，而加速过程的计算量（仅仅做几次四则运算而不需要求函数值）可以忽略不计，可见龙贝格加速过程的效果是极其显著的。

运行示例：

程序源码：

#include <iostream>
#include <cmath>using namespace std;double f(double x)
{return sin(x) / x;
}int main(void)
{double a, b;cout << "请输入积分区间:";cin >> a >> b;double accuracy;cout << "请输入精度:";cin >> accuracy;// 步长double h = b - a;// T1：二分前的梯形法积分值double T1 = h / 2 * (1 + f(b));// T2：二分后的梯形法积分值double T2 = 0;cout << "T" << pow(2, 0) << " = " << T1 << endl;// 对 T1 与 T2 加权平均，得辛普森积分值double S1, S2;// 对 S1 与 S2 加权平均，得柯特斯积分值double C1, C2;// 对 C1 与 C2 加权平均，得龙贝格积分值double R1, R2;// flag 作为循环控制标志性变量int flag = 1;// 记录二分次数int k = 1;while (flag == 1){// 各分点的函数值和double sum = 0;// 分点值double x = a + h / 2;// 在区间上限范围内求各分点的函数值和while (x < b){sum += f(x);x += h;}// 计算梯形序列得下一个二分结果T2 = T1 / 2 + h / 2 * sum;cout << "T" << pow(2, k) << " = " << T2 << endl;// 线性组合外推至辛普生S2 = T2 + 1.0 / 3 * (T2 - T1);cout << "S" << pow(2, k - 1) << " = " << S2 << endl;if (k == 1){// 至少外推 2 次得出 S1，S2k++;h /= 2;T1 = T2;S1 = S2;continue;}else{// 线性组合外推值柯特斯C2 = S2 + 1.0 / 15 * (S2 - S1);cout << "C" << pow(2, k - 2) << " = " << C2 << endl;if (k == 2){// 至少外推 3 次得出C1，C2C1 = C2;k++;h /= 2;T1 = T2;S1 = S2;continue;}else{// 线性组合外推至龙贝格R2 = C2 + 1.0 / 63 * (C2 - C1);cout << "R" << pow(2, k - 3) << " = " << R2 << endl;if (k == 3){// 至少外推 4 次得出 R1, R2R1 = R2;C1 = C2;k++;h /= 2;T1 = T2;S1 = S2;continue;}else if (abs(R2 - R1) >= accuracy){// 精度仍然不符合要求，继续二分步长、继续外推R1 = R2;C1 = C2;k = k + 1;h = h / 2;T1 = T2;S1 = S2;}else{// 精度符合要求，修改 flag 为 0，跳出 while 循环flag = 0;cout << "龙贝格算法求得数值积分结果为:" << R2 << endl;}}}}return 0;
}

[计算机数值积分]龙贝格公式求数值积分相关推荐

龙贝格算法求数值积分的Python程序
Ronberg Integration 分成n等份时辛普森公式的值 from sympy import *def f(t):f = 2000*log(140000/(140000-2100*t))-9 ...
Romberg(龙贝格)求积公式求解数值积分时的注意事项
<数值分析>第5版(李庆扬编著)的第四章课后习题第8-(2)题中,要求使用Romberg(龙贝格)求积公式求解f(x)=xsinx在区间[0,2pi]上的积分,要求误差小于10^(-5). ...
龙贝格公式(一道练习题的解答)
龙贝格公式下标规范: a-d: T=[];% 梯形公式序列 S=[];% 辛普森公式序列 C=[];% 科特斯公式序列 R=[];% 龙贝格公式序列 Y=[]; U=[]; I=[]; S1=0;% ...
数值积分-龙贝格(Romberg)积分
数值积分在工程上是个比较有用的数学工具.在工程上有很多数学问题,看似简单,计算所用的数学公式不算复杂,但是求解起来却很困难,很难获得解析解的公式,这个时候就需要用到数值求解的办法来获取满足工程需要的近 ...
数值积分：龙贝格求积
一.数学原理在变步长的复化梯形计算过程中运用: 就能将粗糙的梯形值Tn逐步加工成精度较高的辛普森值Sn.柯特斯值Cn和龙贝格值Rn.或者说,将收敛缓慢的梯形值序列Tn加工成收敛迅速的龙贝格值序列Rn ...
抛物线求积公式求积分算法matlab,数值计算实验报告---复合求积公式（梯形、抛物线、龙贝格）、导数求值（三点、四点、五点公式）...
----------------------个人作业,如果有后辈的作业习题一致,可以参考学习,一起交流,请勿直接copy ··复合抛物线公式: ··龙贝格公式: 四.实验内容 ------1 实验题目 ...
Romberg(龙贝格)数值积分算法较高效的python实现
1. 原理与公式 2. Python代码实现 # coding=utf-8 import numpy as np import math# 二分法梯形公式 # func:需要积分的函数 # x_min ...
龙贝格数值分析作业c语言,数值分析龙贝格实验报告.doc
数值分析龙贝格实验报告实验三龙贝格方法 [实验类型] 验证性 [实验学时] 2学时 [实验内容] 1.理解龙贝格方法的基本思路 2.用龙贝格方法设计算法,编程求解一个数值积分的问题. [实验前的预 ...
龙贝格算法求解椭球周长
数值积分: 在实际应用中经常应用到计算方法去求解一些不易测量的零件的周长或面积.已知一个椭圆形边框如下图所示,试用龙贝格算法求解这个边框的周长,要求结果精确到6 位有效数字. 3.1 数学原理: 龙贝 ...

[计算机数值积分]龙贝格公式求数值积分

[计算机数值积分]龙贝格公式求数值积分相关推荐

最新文章

热门文章