计算方法实验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}|$小于给定控制精度时停止计算.

代码

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()

结果

部分运算结果

轨迹可视化结果