Romberg积分法数学原理
利用复化梯形求积公式、复化辛普生求积公式、复化柯特斯求积公式的误差估计式计算积分 ∫ a b f ( x ) d x \int_{a}^{b}{f(x)dx} ∫abf(x)dx。记 h = b − a n h = \frac{b - a}{n} h=nb−a, x k = a + k ⋅ h x_{k} = a + k \cdot h xk=a+k⋅h, k = 0 , 1 , ⋯ , n k = 0,1,\cdots,n k=0,1,⋯,n,其计算公式:
T n = 1 2 h ∑ k = 1 n [ f ( x k − 1 ) + f ( x k ) ] T_{n} = \frac{1}{2}h\sum_{k = 1}^{n}\lbrack f(x_{k - 1}) + f(x_{k})\rbrack Tn=21hk=1∑n[f(xk−1)+f(xk)]
T 2 n = 1 2 T n + 1 2 h ∑ k = 1 n f ( x k − 1 2 h ) T_{2n} = \frac{1}{2}T_{n} + \frac{1}{2}h\sum_{k = 1}^{n}{f(x_{k} - \frac{1}{2}h}) T2n=21Tn+21hk=1∑nf(xk−21h)
S n = 1 3 ( 4 T 2 n − T n ) S_{n} = \frac{1}{3}(4T_{2n} - T_{n}) Sn=31(4T2n−Tn)
C n = 1 15 ( 16 S 2 n − S n ) C_{n} = \frac{1}{15}(16S_{2n} - S_{n}) Cn=151(16S2n−Sn)
R n = 1 63 ( 64 C 2 n − C n ) R_{n} = \frac{1}{63}(64C_{2n} - C_{n}) Rn=631(64C2n−Cn)
一般地,利用龙贝格算法计算积分,要输出所谓的 T − T - T−数表
T 1 T 2 S 1 T 4 S 2 C 1 T 8 ⋮ S 4 ⋮ C 2 ⋮ R 1 ⋮ ⋱ \begin{matrix} T_{1} & & & \begin{matrix} & \\ \end{matrix} \\ T_{2} & S_{1} & & \begin{matrix} & \\ \end{matrix} \\ T_{4} & S_{2} & C_{1} & \begin{matrix} & \\ \end{matrix} \\ \begin{matrix} T_{8} \\ \vdots \\ \end{matrix} & \begin{matrix} S_{4} \\ \vdots \\ \end{matrix} & \begin{matrix} C_{2} \\ \vdots \\ \end{matrix} & \begin{matrix} \begin{matrix} R_{1} \\ \vdots \\ \end{matrix} & \begin{matrix} \\ \ddots \\ \end{matrix} \\ \end{matrix} \\ \end{matrix} T1T2T4T8⋮S1S2S4⋮C1C2⋮R1⋮⋱
核心代码
#include <cmath>
#include <cstdio>
#include <iostream>
using namespace std;
#define N 20
int n;
long double a, b, e, T[N][N] = {{0.0}};
long double f(long double x) { return 1.0 / (1.0 + x); }
long double x(int i, long double h) { return a + h * i; }
int main() {
scanf("%llf%llf%llf%d", &a, &b, &e, &n);
int k = 0;
for (; k < n; k++) {
long double h = (b - a) / pow(2, k), sum = 0.0;
for (int i = 1; i <= pow(2, k) - 1; i++) sum += f(x(i, h));
T[k][0] = 0.5 * h * (f(a) + 2 * sum + f(b));
for (int m = 1; m <= k; m++)
T[k][m] =
(pow(4, m) * T[k][m - 1] - T[k - 1][m - 1]) / (pow(4, m) - 1);
for (int m = 0; m <= k; m++) printf("%.7Lf\t", T[k][m]);
cout << endl;
if (k > 0 && fabs(T[k][0] - T[k][k]) <= e) break;
}
return 0;
}