用C# 写的龙贝格（Romberg）积分法

很久之前了，群上有人问有没有用C#写的积分函数，因为我自己以前做过一个积分计算器，就跟他说搜索一下Romberg积分法吧，但是后来他说网上写的都是用C++写的方法，没有C#的，问我有没有，但是我之前是在家里做的，学校里没有副本，可惜。

今天翻一下邮箱，发现居然是我忘记贴标签，导致当时没有找到，感觉有点对不起那位兄弟，在此就拿出来分享一下吧。

网上翻过很多代码都是照搬算法说明来迭代的，我却把那个公式代进去，化简到只剩下四项，将空间要求降低了，如果还有什么能改进的地方请各位多多指教！以下是代码片断：

Code
/// <summary>
/// Romberg 龙贝格算法
/// </summary>
/// <param name="f">被积表达式函数</param>
/// <param name="a">积分下限</param>
/// <param name="b">积分上限</param>
/// <param name="e">精度</param>
/// <param name="maxSteps">最大迭代次数</param>
/// <returns></returns>
public static double Romberg( Func<double, double> f, double a, double b, double e, int maxSteps ) {
    double M = 2835.0;
    double A = 217.0;
    double B = 2048.0;
    double C = 352.0;
    double D = 218.0;

    double T1n = (b - a) * (f(b) + f(a)) / 2;
    double H4n = H(f, a, b, 4);
    double H2n = H(f, a, b, 2);
    double H1n = H(f, a, b, 1);

    double R1 = (A * T1n + B * H4n + C * H2n + D * H1n) / M;
    double R2 = 0.0;

    for (int steps = 0; steps < maxSteps; ++steps) {
        T1n = (T1n + H1n) / 2;
        H1n = H2n;
        H2n = H4n;
        H4n = H(f, a, b, 8 << steps);
        R2 = (A * T1n + B * H4n + C * H2n + D * H1n) / M;

        if (Math.Abs(R2 - R1) <= e) {
            return R2;
        }
        else {
            R1 = R2;
        }
    }
    return R2;
}

private static double H( Func<double, double> f, double a, double b, int n ) {
    double h = (b - a) / n;
    double sum = 0.0;
    for (int k = 0; k < n; ++k) {
        sum += f(a + (k + 0.5) * h);
    }
    return h * sum;
}