如何使用自适应正交例程整合向量的内容

Question

我有一个函数，它要求返回类型是一个容器。问题是我需要尽可能有效地集成容器的内容，并希望使用自适应 Gauss-Kronrod 集成或同样有效（或更好）的方法。

我希望使用 GNU Scientific Library 正交例程 qags，但它返回的结果是 double 类型。就目前情况而言，我认为我可能不得不重写 GSL 正交例程的部分以将返回类型转换为 std 向量，但这将是一个相当冗长且可能容易出错的弯路。我希望有人可以推荐更好的解决方案！

下面是我试图整合的那种函数的例子，目前使用基本的梯形规则，但表明我更喜欢在哪里实现 GSL 例程 gaq。虽然它比我的实际问题简单得多，但它表明放入向量中的每个元素都是根据先前的结果计算的，因此需要容器。

#include <iostream>
#include <vector>
#include <gsl/gsl_integration.h>
using namespace std;

vector<double> f(double E, int N) {
    vector<double> result;
    result.reserve(N);
    double x = E;
    for (int it=0; it < N; ++it){
        result.push_back(x);
        x *= x;
    }
    return result;
}

vector<double> f_qag(double E, void * params) {
    int N = *(int *) params;
    vector<double> result;
    result.reserve(N);
    double x = E;
    for (int it=0; it < N; ++it){
        result.push_back(x);
        x *= x;
    }
    return result;
}


int main(){
    vector<double> result;
    vector<double> integrate;
    int N = 20;
    result.reserve(N);
    integrate.reserve(N);
    for (int i = 0; i < N; i++)
        result[i] = 0.;
    double end = 1.0;
    double start = -1.0;

    // I would like to use qag here
    /* gsl_integration_workspace * w = gsl_integration_workspace_alloc (1000); */
    /* vector<double> error; */
    /* gsl_function F; */
    /* F.function = &f_qag; */
    /* F.params = &N; */
    /* gsl_integration_qag (&F, start, end, 0, 1e-7, 1000, 1, w, &result, &error); */

    // instead of resorting to the trapezoidal rule here
    double E;
    const int n = 1000;
    double factor = (end - start)/(n*1.);
    for (int k=0; k<n+1; k++) {
        E = start + k*factor;
        integrate = f(E, N);
        for (int i = 0; i < N; i++){
            if ((k==0)||(k==n))
                result[i] += 0.5*factor*integrate[i];
            else 
                result[i] += factor*integrate[i];
        }
    }   

    for (int i = 0; i < N; i++)
        cout<<i<<" "<<result[i]<<endl;

    return 0;
}

我打算检查与伪结果 conv 的收敛性；

double conv = 0.;
for (int i = 0; i < N; i++)
    conv += abs(integrate[i]);