基于离散点的快速C++数值积分库求助

• 需对三维空间执行大量积分运算
• 已获取空间的均匀离散网格（如各方向30个点），仅已知网格点的待积分函数值，无解析函数形式 f(x,y,z)
• 原使用MATLAB的trapz(...)实现，现计划迁移至C++以减少单积分开销、加速代码
• GNU科学库（GSL）因要求函数句柄及解析形式，无法满足需求

1. Eigen #

轻量级线性代数库，可手动实现多维度梯形积分（类似trapz），速度快、无额外依赖，适合自定义离散积分逻辑：

#include <Eigen/Dense>

double trapz3D(const Eigen::ArrayXd& x, const Eigen::ArrayXd& y, const Eigen::ArrayXd& z, const Eigen::Array<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::Dynamic>& f) {
    double dx = x(1) - x(0);
    double dy = y(1) - y(0);
    double dz = z(1) - z(0);
    
    // 三维梯形积分：先对z积分，再y，再x
    Eigen::ArrayXXd f_z = f.trapz(z);
    Eigen::ArrayXd f_zy = f_z.trapz(y);
    double integral = f_zy.trapz(x);
    
    return integral * dx * dy * dz;
}

注：Eigen 3.3及以上版本的Array类内置trapz方法，直接支持离散数据的单维度积分，嵌套调用即可实现三维积分，性能接近手写优化代码。

2. Boost.Math #

提供成熟的数值积分工具，支持离散数据的梯形积分，兼容多种容器类型，也可适配后续非均匀网格需求：

#include <boost/math/quadrature/trapezoidal.hpp>
#include <vector>

using namespace boost::math::quadrature;

double trapz3D(const std::vector<double>& x, const std::vector<double>& y, const std::vector<double>& z, const std::vector<std::vector<std::vector<double>>>& f) {
    double dx = x[1] - x[0];
    double dy = y[1] - y[0];
    double dz = z[1] - z[0];
    
    // 逐层积分
    std::vector<std::vector<double>> f_z(f.size(), std::vector<double>(f[0].size()));
    for (size_t i = 0; i < f.size(); ++i) {
        for (size_t j = 0; j < f[0].size(); ++j) {
            f_z[i][j] = trapezoidal(f[i][j].begin(), f[i][j].end(), dz);
        }
    }
    
    std::vector<double> f_zy(f.size());
    for (size_t i = 0; i < f.size(); ++i) {
        f_zy[i] = trapezoidal(f_z[i].begin(), f_z[i].end(), dy);
    }
    
    double integral = trapezoidal(f_zy.begin(), f_zy.end(), dx);
    return integral;
}

3. Fastor #

专为张量运算优化的C++库，通过编译期优化生成高效代码，适合固定网格尺寸下的极致性能需求：

#include <Fastor/Fastor.h>

using namespace Fastor;

template<size_t I, size_t J, size_t K>
double trapz3D(const Tensor<double,I>& x, const Tensor<double,J>& y, const Tensor<double,K>& z, const Tensor<double,I,J,K>& f) {
    double dx = x(1) - x(0);
    double dy = y(1) - y(0);
    double dz = z(1) - z(0);
    
    // 利用切片和求和实现梯形积分逻辑
    Tensor<double,I,J> f_z = 0.5 * dz * (f.slice(K,0) + f.slice(K,K-1) + 2.0 * sum(f.slice(K,1,K-1), Axis<K>()));
    Tensor<double,I> f_zy = 0.5 * dy * (f_z.slice(J,0) + f_z.slice(J,J-1) + 2.0 * sum(f_z.slice(J,1,J-1), Axis<J>()));
    double integral = 0.5 * dx * (f_zy.slice(I,0) + f_zy.slice(I,I-1) + 2.0 * sum(f_zy.slice(I,1,I-1), Axis<I>()));
    
    return integral;
}

• 轻量无依赖需求：选Eigen，代码简洁且性能足够
• 已有Boost生态或需兼容非均匀网格：选Boost.Math
• 固定网格下追求极致性能：选Fastor

内容的提问来源于stack exchange，提问作者blockchain187