计算方法实验9：Romberg积分求解速度、位移
计算方法实验9：Romberg积分求解速度、位移
本文介绍了使用Romberg积分方法求解质点运动轨迹的问题。通过详细描述任务要求、算法原理、代码实现以及结果展示，帮助读者理解数值积分在实际问题中的应用。
任务
- 输出质点的轨迹$(x(t), y(t))$, $t \in {0.1, 0.2, 0.3, ..., 10}$，并在二维平面中画出该轨迹。
- 请比较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(f)-S_{n}(f)=-\frac{f^{(4)}(\xi)}{180}h^{4}(b-a)\approx dh^{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(f)\approx S_{2n}(f)+\frac{1}{15}\left(S_{2n}(f)-S_{n}(f)\right)=C_{n}(f)$$
复化 Simpson 公式组成复化 Cotes 公式，其截断误差是$O(h^6)$。同理对 Cotes公式进行线性组合：
$$I(f)-C_{2n}(f)=e\left(\frac{h}{2}\right)^{6}$$
$$I(f)-C_{n}(f)=eh^{6}$$
得到具有 7 次代数精度和截断误差是$O(h^8)$的 Romberg 公式：
$$R_{n}(f)=C_{2n}(f)+\frac{1}{63}\left(C_{2n}(f)-C_{n}(f)\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}|$小于给定控制精度时停止计算。
代码
C++实现Romberg积分运算
#include<bits/stdc++.h>
using namespace std;
int satisfiedCount;
int M;
long double eps = 1e-6, proportion;
vector<long double> Vx(100), Vy(100);
long double ax(long double t);
long double ay(long double t);
long double vx(long double t);
long double vy(long double t);
long double romberg(function<long double(long double)> f, long double a, long double b, long double eps, int M, bool isX);// Perform the Romberg integration
int main()
{
for (M = 4; M <= 20; M += 4)
{
satisfiedCount = 0;
cout << "M = " << M << endl;
ofstream outFile("trajectory.txt", ios::app);
for (long double t = 0.1; t < 10.1; t += 0.1)
{
long double v_x = vx(t);
long double v_y = vy(t);
long double x = romberg(vx, 0.0, t, eps, M, 1);
long double y = romberg(vy, 0.0, t, eps, M, 1);
cout << fixed << setprecision(1) << "At t = " << t << ", vx = " << setprecision(6) << v_x << ", vy = " << setprecision(6) << v_y << ", " << "(x, y) = (" << setprecision(6) << x << ", " << setprecision(6) << y << ")" << endl;
if (M == 8)
outFile << fixed << setprecision(6) << x << " " << y << "\n";
}
long double proportion = (long double)satisfiedCount / 200;
cout << "At M = " << M << ", proportion of times the error requirement of (x,y) 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);
}
long double vx(long double t)
{
return romberg(ax, 0.0, t, eps, M, 0);
}
long double vy(long double t)
{
return romberg(ay, 0.0, t, eps, M, 0);
}
// Perform the Romberg integration
long double romberg(function<long double(long double)> f, long double a, long double b, long double eps, int M, bool isX) {
long double h[M+1], r[M+1][M+1];
h[1] = b - a;
r[1][1] = 0.5 * h[1] * (f(a) + f(b));
for (int k = 2; k <= M; k++)
{
h[k] = 0.5 * h[k-1];
long double sum = 0;
for (int i = 1; i <= pow(2, k-2); i++)
sum += f(a + (2*i-1) * h[k]);
r[k][1] = 0.5 * (r[k-1][1] + h[k-1] * sum);
for (int j = 2; j <= k; j++)
r[k][j] = r[k][j-1] + (r[k][j-1] - r[k-1][j-1]) / (pow(4, j-1) - 1);
if (k > 2 && fabs(r[k][k] - r[k-1][k-1]) < eps)
{
if(isX)
satisfiedCount++;
return r[k][k];
}
}
return r[M][M];
}
Python可视化运动轨迹
import matplotlib.pyplot as plt
with 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()