计算方法实验9：Romberg积分求解速度、位移

任务

1. 输出质点的轨迹 ( x ( t ) , y ( t ) ) , t ∈ { 0.1 , 0.2 , 0.3 , . . . , 10 } (x(t), y(t)), t\in \{0.1, 0.2, 0.3, ..., 10\} (x(t),y(t)),t∈{0.1,0.2,0.3,...,10}，并在二维平面中画出该轨迹.
2. 请比较M分别取4, 8, 12, 16, 20 时，Romberg积分达到要求精度的比例（达到误差要求的次数/调用总次数），分析该比例随M 的变化。

算法

现在要用数值方法求${\int }_a^{b}f(x)dx$，设$h=\frac{b-a}{n}$，已知：

复化梯形积分$T_n\left(f\right)=h\left[\frac{1}{2}f\left(a\right)+{\sum }_{i=1}^{n-1}f\left(a+ih\right)+\frac{1}{2}f\left(b\right)\right]$、

复化Simpson积分$S_n\left(f\right)=\frac{h}{3}\left[f\left(a\right)+4{\sum }_{i=0}^{m-1}f\left(x_{2i+1}\right)+2{\sum }_{i=1}^{m-1}f\left(x_{2i}\right)+f\left(b\right)\right]$.

将 $(T_n(f)-T_{2n}(f))$ 作 为 $T_{2n}(f)$的修正值补充到$I(f)$，即
$$I(f)\approx T_{2n}(f)+\frac{1}{3}\left(T_{2n}\left(f\right)-T_n\left(f\right)\right)=\frac{4}{3}{T}_{2n}-\frac{1}{3}{T}_n=S_n$$

其结果是将复化梯形求积公式组合成复化 Simpson 求积公式， 截断误差由$O(h^2)$提高到$O(h^4)$，这种手段称为外推算法，该算法在不增加计算量的前提下提高了误差的精度.不妨对$S_{2n}(f),S_n(f)$再作一次线性组合：

$$I\left(f\right)-S_n\left(f\right)=-\frac{f^{\left(4\right)}\left(\xi \right)}{180}{h}^4\left(b-a\right)\approx d{h}^4$$
$$I(f)-S_{2n}(f)=-\frac{f^{(4)}(\eta )}{180}{\left(\frac{h}{2}\right)}^4(b-a)\approx d{\left(\frac{h}{2}\right)}^4$$

$$I\left(f\right)\approx S_{2n}\left(f\right)+\frac{1}{15}\left(S_{2n}\left(f\right)-S_n\left(f\right)\right)=C_n\left(f\right)$$
复化 Simpson 公式组成复化 Cotes 公式，其截断误差是$O(h^6).$同理对 Cotes公式进行线性组合：
$$\begin{gathered} I\left(f\right)-C_{2n}\left(f\right)=e{\left(\frac{h}{2}\right)}^6 \\ I\left(f\right)-C_n\left(f\right)=e{h}^6 \end{gathered}$$
得到具有 7 次代数精度和截断误差是$O(h^8)$的 Romberg 公式：
$$R_n\left(f\right)=C_{2n}\left(f\right)+\frac{1}{63}\left(C_{2n}\left(f\right)-C_n\left(f\right)\right)$$

为了便于在计算机上实现 Romberg 算法，将$T_n,S_n,C_n,R_n,\cdots$统一用$R_{k,j}$ 表示，列标$j=1,2,3,4$分别表示梯形、Simpson、Cotes 、Romberg积分，行标$k$表示步长$h_k=\frac{h}{2^{k-1}}$，得到Romberg 计算公式：
$$R_{k,j}=R_{k,j-1}+\frac{R_{k,j-1}-R_{k-1,j-1}}{4^{j-1}-1},k=2,3,\cdots$$
对每一个$k,j$从 2 做到$k$,一直做到$|R_k,k-R_{k-1,k-1}|$小于给定控制精度时停止计算.
注：下面代码中数组下标从0开始.

代码

C++实现Romberg积分运算

#include<bits/stdc++.h>
using namespace std;int satisfiedCount;long double ax(long double t);
long double ay(long double t);
long double romberg(function<long double(long double)> f, long double a, long double b, long double eps, int maxIter, bool isX);// Perform the Romberg integrationint main() 
{long double eps = 1e-6, proportion;int maxIter;satisfiedCount = 0;maxIter = 4;cout << "maxIter = " << maxIter << endl;for (long double t = 0.1; t <= 10; t += 0.1) { long double vx = romberg(ax, 0, t, eps, maxIter, 0);long double vy = romberg(ay, 0, t, eps, maxIter, 0);long double x = romberg([&](long double t) { return vx; }, 0, t, eps, maxIter, 1);long double y = romberg([&](long double t) { return vy; }, 0, t, eps, maxIter, 0);cout << fixed << setprecision(1) << "At t = " << t << ", vx = " << setprecision(6) << vx << ", vy = " << setprecision(6) << vy << ", " << "(x, y) = (" << setprecision(6) << x << ", " << setprecision(6) << y << ")" << endl;}proportion = (long double)satisfiedCount / 100;cout << "At maxIter = " << maxIter << ", proportion of times the error requirement of x was satisfied: " << proportion << endl;satisfiedCount = 0;maxIter = 8;cout << "maxIter = " << maxIter << endl;ofstream outFile("trajectory.txt");for (long double t = 0.1; t <= 10; t += 0.1) {long double vx = romberg(ax, 0, t, eps, maxIter, 0);long double vy = romberg(ay, 0, t, eps, maxIter, 0);long double x = romberg([&](long double t) { return vx; }, 0, t, eps, maxIter, 1);long double y = romberg([&](long double t) { return vy; }, 0, t, eps, maxIter, 0);cout << fixed << setprecision(1) << "At t = " << t << ", vx = " << setprecision(6) << vx << ", vy = " << setprecision(6) << vy << ", " << "(x, y) = (" << setprecision(6) << x << ", " << setprecision(6) << y << ")" << endl;outFile << fixed << setprecision(6) << x << " " << y << "\n";//把坐标写入文件，方便画轨迹}proportion = (long double)satisfiedCount / 100;cout << "At maxIter = " << maxIter << ", proportion of times the error requirement of x was satisfied: " << proportion << endl;satisfiedCount = 0;maxIter = 12;cout << "maxIter = " << maxIter << endl;for (long double t = 0.1; t <= 10; t += 0.1) {long double vx = romberg(ax, 0, t, eps, maxIter, 0);long double vy = romberg(ay, 0, t, eps, maxIter, 0);long double x = romberg([&](long double t) { return vx; }, 0, t, eps, maxIter, 1);long double y = romberg([&](long double t) { return vy; }, 0, t, eps, maxIter, 0);cout << fixed << setprecision(1) << "At t = " << t << ", vx = " << setprecision(6) << vx << ", vy = " << setprecision(6) << vy << ", " << "(x, y) = (" << setprecision(6) << x << ", " << setprecision(6) << y << ")" << endl;}proportion = (long double)satisfiedCount / 100;cout << "At maxIter = " << maxIter << ", proportion of times the error requirement of x was satisfied: " << proportion << endl;satisfiedCount = 0;maxIter = 16;cout << "maxIter = " << maxIter << endl;for (long double t = 0.1; t <= 10; t += 0.1) {long double vx = romberg(ax, 0, t, eps, maxIter, 0);long double vy = romberg(ay, 0, t, eps, maxIter, 0);long double x = romberg([&](long double t) { return vx; }, 0, t, eps, maxIter, 1);long double y = romberg([&](long double t) { return vy; }, 0, t, eps, maxIter, 0);cout << fixed << setprecision(1) << "At t = " << t << ", vx = " << setprecision(6) << vx << ", vy = " << setprecision(6) << vy << ", " << "(x, y) = (" << setprecision(6) << x << ", " << setprecision(6) << y << ")" << endl;}proportion = (long double)satisfiedCount / 100;cout << "At maxIter = " << maxIter << ", proportion of times the error requirement of x was satisfied: " << proportion << endl;satisfiedCount = 0;maxIter = 20;cout << "maxIter = " << maxIter << endl;for (long double t = 0.1; t < 10.1; t += 0.1) {long double vx = romberg(ax, 0, t, eps, maxIter, 0);long double vy = romberg(ay, 0, t, eps, maxIter, 0);long double x = romberg([&](long double t) { return vx; }, 0, t, eps, maxIter, 1);long double y = romberg([&](long double t) { return vy; }, 0, t, eps, maxIter, 0);cout << fixed << setprecision(1) << "At t = " << t << ", vx = " << setprecision(6) << vx << ", vy = " << setprecision(6) << vy << ", " << "(x, y) = (" << setprecision(6) << x << ", " << setprecision(6) << y << ")" << endl;}proportion = (long double)satisfiedCount / 100;cout << "At maxIter = " << maxIter << ", proportion of times the error requirement of x was satisfied: " << proportion << endl;return 0;
}long double ax(long double t) 
{return sin(t) / (1 + sqrt(t));
}long double ay(long double t) 
{return log(t + 1) / (t + 1);
}// Perform the Romberg integration
long double romberg(function<long double(long double)> f, long double a, long double b, long double eps, int maxIter, bool isX) {long double h[maxIter], r[maxIter][maxIter];h[0] = b - a;r[0][0] = 0.5 * h[0] * (f(a) + f(b));for (int i = 1; i < maxIter; i++) {h[i] = 0.5 * h[i-1];long double sum = 0;for (int k = 0; k < pow(2, i-1); k++)sum += f(a + (2*k+1) * h[i]);r[i][0] = 0.5 * r[i-1][0] + h[i] * sum;for (int j = 1; j <= i; j++)r[i][j] = r[i][j-1] + (r[i][j-1] - r[i-1][j-1]) / (pow(4, j) - 1);if (i > 1 && fabs(r[i][i] - r[i-1][i-1]) < eps){if(isX)satisfiedCount++;return r[i][i];}}return r[maxIter-1][maxIter-1];
}

python可视化运动轨迹

import matplotlib.pyplot as pltwith open('trajectory.txt', 'r') as file:lines = file.readlines()x, y = zip(*[(float(line.split()[0]), float(line.split()[1])) for line in lines])plt.plot(x, y, 'o-')
plt.xlabel('X')
plt.ylabel('Y')
plt.title('Plot of points with smooth curve')
plt.show()

结果

部分运算结果

轨迹可视化结果