Как интегрировать содержимое вектора с помощью адаптивной квадратурной процедуры - PullRequest
Новая
       36

Как интегрировать содержимое вектора с помощью адаптивной квадратурной процедуры

0 голосов
/ 26 апреля 2019

У меня есть функция, которая требует, чтобы тип возвращаемого значения был контейнером. Проблема в том, что мне нужно максимально эффективно интегрировать содержимое контейнера, и я надеялся использовать адаптивную интеграцию Гаусса-Кронрода или что-то такое же эффективное (или лучшее).

Я надеялся использовать квадратурную процедуру qags научной библиотеки GNU, но она возвращает результат с типом 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]);
...