c++ - 在C++中实现一维gsl集成,参数作为集成边界
问题描述
我想使用 gsl 在 C++ 中执行以下一维积分,
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++ 的新手,所以提前为我的问题的简单性道歉。
解决方案
如果积分下限(“x/4”表达式)中表示的 x 与参与被积函数公式的 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 . . .
推荐阅读
- python - Python - 从列表中搜索数据框中的字符串
- queue - 在双端队列中,我们可以在两端插入和删除。它是否违反队列属性?
- javascript - 无法理解 javascript 函数的输出
- search - 如果使用欧几里得距离启发式的 A* 搜索允许对角线移动,它仍然是最优的吗?
- java - 注释和方法调用之间的 Mockito 区别
- python - 使用 tf.data 实现 tensorflow 输入管道时出错
- excel - 如何在 Excel 中创建部分字符串搜索功能,查看整个数组
- c - 如何化解这个二元炸弹阶段 4
- linux - 使用 NVidia GPU 的屏幕 EGL 显示
- javascript - 无法在网格 Sencha Extjs 中选择字段