是否可以优化此集成代码使其运行得更快？

计算科学 C++ 表现
2021-12-05 06:41:23 
double trap(double func(double), double b, double a, double N) {
  double j;
  double s;
  double h = (b-a)/(N-1.0); //Width of trapezia

  double func1 = func(a);
  double func2;

  for (s=0,j=a;j<b;j+=h){
    func2 = func(j+h);
    s = s + 0.5*(func1+func2)*h;
    func1 = func2;
  }

  return s;
}

以上是我的 C++ 代码，用于使用func()梯形在限制之间进行一维数值积分（使用扩展梯形规则）。[a,b]N−1

我实际上是在做一个 3D 集成，这个代码是递归调用的。我使用给我带来了不错的结果。N=50

除了进一步减少之外，是否有人能够建议如何优化上面的代码以使其运行得更快？或者，甚至可以建议一种更快的集成方法？N

4个回答

函数的评估很可能是该计算中最耗时的部分。如果是这种情况，那么您应该专注于提高 func() 的速度，而不是试图加快集成例程本身。

根据 func() 的属性，您还可能通过使用更复杂的积分公式以更少的函数评估获得更精确的积分评估。

可能的？是的。有用？不。我将在此处列出的优化不太可能在运行时产生超过百分之几的差异。一个好的编译器可能已经为你做了这些。

无论如何，看看你的内部循环：

    for (s=0,j=a;j<b;j+=h){
        func2 = func(j+h);
        s = s + 0.5*(func1+func2)*h;
        func1 = func2;
    }

在每次循环迭代中，您执行三个可以带入外部j + h的数学运算：加法、乘法0.5和乘法h。第一个您可以通过在 处启动您的迭代器变量来解决，a + h其他的可以通过分解乘法来解决：

    for (s=0, j=a+h; j<=b; j+=h){
        func2 = func(j);
        s += func1+func2;
        func1 = func2;
    }
    s *= 0.5 * h;

虽然我会指出，通过这样做，由于浮点舍入错误，可能会错过循环的最后一次迭代。（这也是您原始实现中的一个问题。）要解决这个问题，请使用unsigned intorsize_t计数器：

    size_t n;
    for (s=0, n=0, j=a+h; n<N; n++, j+=h){
        func2 = func(j);
        s += func1+func2;
        func1 = func2;
    }
    s *= 0.5 * h;

正如布赖恩的回答所说，您最好将时间花在优化函数的评估上func。如果这种方法的准确性足够，我怀疑你会发现同样的方法更快N。（尽管您可以运行一些测试来查看例如 Runge-Kutta 是否可以让您降低N到足以使整体集成花费更少时间而不牺牲准确性的程度。）

在数学上，您的表达式相当于：

I=h(12f1+f2+f3+...+fn−1+12fn)+O((b−a)3f″n2)

所以你可以实现它。如前所述，时间可能主要由函数评估决定，因此要获得相同的准确度，您可以使用需要较少函数评估的更好的集成方法。

在现代，高斯求积不仅仅是一个玩具。仅当您需要很少的评估时才有用。如果你想要一些容易实现的东西，你可以使用辛普森的规则，但如果没有充分的理由，我不会比 order1/N3

如果函数的曲率变化很大，您可以使用自适应步长例程，当函数平坦时选择较大的步长，曲率较高时选择较小的更精确的步长。

我建议进行一些更改以改进计算：

  • 为了性能和准确性，请使用std::fma()，它执行融合乘加。
  • 为了性能，推迟将每个梯形的面积乘以 0.5——你可以在最后做一次。
  • 避免重复添加h，这可能会累积舍入误差。

此外，为了清楚起见，我会进行一些更改：

  • 给函数一个更具描述性的名称。
  • 交换函数签名中a和的顺序。b
  • 重命名N→ n、h→ dx、j→ x2、s→ accumulator。
  • 更改n为int.
  • 在更严格的范围内声明变量。
#include <cmath>

double trapezoidal_integration(double func(double), double a, double b, int n) {
    double dx = (b - a) / (n - 1);   // Width of trapezoids

    double func_x1 = func(a);
    double accumulator = 0;

    for (int i = 1; i <= n; i++) {
        double x2 = a + i * dx;      // Avoid repeated floating-point addition
        double func_x2 = func(x2);
        accumulator = std::fma(func_x1 + func_x2, dx, accumulator); // Fused multiply-add
        func_x1 = func_x2;
    }

    return 0.5 * accumulator;
}
其它你可能感兴趣的问题
上一篇使用 Armadillo 对矩阵乘法进行超级 C++ 优化 下一篇哪种 Runge-Kutta 方法更准确：Dormand-Prince 或 Cash-Karp？