Реализация одномерной интеграции gsl в C++ с параметрами в качестве границы интеграции

Я хотел бы выполнить следующую одномерную интеграцию в С++ с использованием gsl,

I(x) = int_{x/4}^1 dy y^{3/2+a} (1-y)^{1/2} exp(1.13*sqrt(log(4y/x))),, где a и x рассматриваются как константы интегрирования, и я хотел бы иметь возможность выполнить интегрирование для x, выбранных в интервале 0.000001<x<0.001 для фиксированного a, предоставленного пользователем.

Вот моя попытка:

 #include <iostream>
 #include <iomanip>
 #include <fstream>
 #include <vector>
 #include <string>
 #include <cmath>
 #include <gsl/gsl_integration.h>
 #include <stdio.h>
 #include <math.h>

double integrand(double y, void * params) {
double a = *(double *) params;
double x = *(double *) params; 

double intg = pow(y,3e0/2e0+a)*pow(1-y,1e0/2e0)*exp(1.13*sqrt(log(4*y/x)));

return intg;

}

double integral(double a) {

gsl_integration_workspace * w
= gsl_integration_workspace_alloc (1000);
gsl_function F;
F.function = &integrand;                 
F.params = &a;
double result, error;

for(int i = 0; i<100; i++) {
    double x = (0.001-0.000001)*i/100 + 0.000001;


gsl_integration_qags (&F, x/4, 1e0, 0, 1e-7, 1000,
                      w, &result, &error);

}
 gsl_integration_workspace_free (w); // Free memory

return result;
}




int main(){
std::cout << "x" << x << "result "<< integral(-0.046)<<std::endl;

 }

Проблема, с которой я сталкиваюсь, заключается в том, как передать значение x, заданное циклом for, в подынтегральное выражение? На данный момент код возвращает nan, потому что я не передаю x в подынтегральное выражение, а только в зависимость от x в нижней границе интегрирования. Все еще новичок в C++, поэтому заранее извиняюсь за возможную простоту моего вопроса.


person CAF    schedule 21.05.2020    source источник


Ответы (1)


Если x, указанный в нижнем пределе интегрирования (выражение x/4), совпадает с x, участвующим в формуле подынтегрального выражения, то решение будет следующим:

#include <iostream>
#include <iomanip>
#include <fstream>
#include <vector>
#include <string>
#include <cmath>
#include <gsl/gsl_integration.h>
#include <stdio.h>
#include <math.h>

//create a structure for the parameters of the integral function: 

struct params 
{
    double a;
    double x;
};

//create an integration function for the x variable: 
//let all INTERNAL parameters have a prefix "_" in the name:   

double integrand(double y, void * params_from_above)
{
    //declare our structure for parameters and be sure to do type casting, 
    //and initialize the parameter values passed from above:
    params& _inner_params = *(params*)params_from_above;
    
    //we calculate the value of the integrand using the passed parameter values:
    double _intg = pow(y, 1.5 + _inner_params.a) * pow(1.0 - y, 0.5) * exp(1.13 * sqrt(log(4.0 * y/_inner_params.x)));

    return _intg;
}


double integral(double a)
{
    double result, error;
    
    //declare a structure for storing and passing "a" and "x" parameters:
    params custom_params;
    //create integration workspace:
    gsl_integration_workspace * w = gsl_integration_workspace_alloc (1000);
    
    //declare and initialize gsl callback function for integrand:
    gsl_function F;
    F.function = &integrand;
    
    //declare and initialize gsl params structure:
    //(so far only for the value of "a" parameter)
    F.params = &custom_params;    
    custom_params.a = a;

    //start the main cycle of calculations:    
    for(int i = 0; i < 100; i++)
    {
        //compute x parameter value:
        double _x = (0.001-0.000001)*i/100 + 0.000001;
        //save the value of the X parameter for sending to the integrand:
        custom_params.x = _x;
        
        //compute the integral for current "x"/4 parameter value:
        gsl_integration_qags (&F, _x/4.0, 1.0, 0, 1e-7, 1000, w, &result, &error);
        
        //output integral value:
        printf("res = %f\n", result);
        //delay: press any key for continue: 
        system("pause");

    }
    gsl_integration_workspace_free (w); // Free memory

    return result;
}

int main()
{
    //here you are mistaken, since only the last value is displayed...
    //but I left your expression for the sake of calling an integral function:
    std::cout << "result "<< integral(-0.046) << std::endl; 
} 

И я получаю следующий вывод для первых пяти значений параметра x:

C:\Users\...\CLionProjects\gsl_test\cmake-build-debug\gsl_test.exe
res = 15.226887
Press any key to continue . . .

res = 10.523077
Press any key to continue . . .

res = 9.467232
Press any key to continue . . .

res = 8.870463
Press any key to continue . . .

res = 8.459477
Press any key to continue . . .
person DJNZ    schedule 30.06.2020
comment
большое спасибо, так и есть. Не могли бы вы объяснить, что означает фрагмент кода «params& _inner_params = (params)params_from_above;’? Я понимаю, что выполняется некоторое приведение, потому что внутренняя переменная является указателем на тип void, а элементы params имеют тип double? Но я хотел бы понять этот фрагмент по частям. Большое спасибо. - person CAF; 01.07.2020
comment
params& _inner_params = *(params*)params_from_above; слева: что означает, что мы создаем ссылку на тип params (ссылка на тип структуры). Справа: мы конвертируем типы явно: преобразуем тип void* в тип params* (тип указателей в тип указателей!) и берем значения из указателя! Это очень простое выражение, читайте его справа налево: возьмите указатель void, преобразуйте в указатель params, затем возьмите значение из этого указателя и инициализируйте ссылку этим значением. - person DJNZ; 02.07.2020
comment
Спасибо. Итак, чтобы проверить мое понимание: params_from_above - это указатель типа void, поэтому *(params*)params_from_above приводит params_from_above к типу params. Затем внутренняя переменная _inner_params инициализируется ЗНАЧЕНИЕМ, на которое указывает этот указатель, с помощью функции &. Это правильное значение? - person CAF; 14.07.2020
comment
Да, почти правильно, но проще :) _inner_params -- это просто референция, то есть псевдоним, псевдоним для громоздкой конструкции *(params*)params_from_above. Вместо ссылки используйте *(params*)params_from_above, потому что ссылка -- проще, короче и безопаснее. Еще раз — объявляем _inner_params как переменную типа структуры и тут же инициализируем, беря ЗНАЧЕНИЕ из указателя (ЗНАЧЕНИЕ берем с помощью крайнего левого выражения * в *(params*)params_from_above). P.S.: если моя версия кода вам помогла, отметьте ее как решение :) - person DJNZ; 15.07.2020
comment
большое спасибо, это действительно имеет больше смысла. +1 и принято - person CAF; 15.07.2020