GSL----蒙特卡洛多重积分

最新推荐文章于 2024-03-31 16:21:23 发布
最新推荐文章于 2024-03-31 16:21:23 发布
阅读量790 收藏
点赞数
文章标签： matlab
原文链接：http://www.cnblogs.com/JustHaveFun-SAN/archive/2011/03/23/san.html
版权

GSL-蒙特卡洛积分  TR:SAN E:VISUALSAN@YAHOO.CN  2011.3.23

---------------------------------------------------------------------------

Gsl中包含有计算N重积分的蒙特卡洛方法，不过只能计算积分上下限确定的多重积分，对于上下限带有参数的积分无能为力，表达式如下：

三种蒙特卡洛积分实现方法分别定义在头文件：

#include <gsl/gsl_monte_miser.h>

#include <gsl/gsl_monte_plain.h>

#include <gsl/gsl_monte_vegas.h>

每组方法的计算思路是一样的：生成控制变量、初始化控制变量、计算积分、释放控制变量。每种方法都将用到随机数产生器。其中需要用户自定义被积分函数：

double (*f)(double * x_array, size_t dim, void * params);

其中x是积分变量数组，dim为积分变量个数，parms为额外传递的参数

对于二次函数f(x,y)=ax2+bxy+cy2，当a=3,b=2,c=1时，下面将定义目标函数：

struct my_par
{
 double a,b,c;
};
 double my_f( double * x_array, size_t dim, void * params )
{
my_par * mp = (my_par * ) params ;
 if (dim != 2 )
{
fprintf(stderr, " dim!=2 " );
abort();
}
 return x_array[ 0 ] * x_array[ 0 ] * mp -> a +
x_array[ 0 ] * x_array[ 1 ] * mp -> b +
x_array[ 1 ] * x_array[ 1 ] * mp -> c;
}

 void test_monte()
{
my_par par = { 3 , 2 , 1 };
gsl_monte_function f;
f.dim = 2 ;
f.f = my_f;
f. params =& par;

}

下面介绍用法：

  1. plainmonte carlo

     包含头文件：#include <gsl/gsl_monte_plain.h>

     使用步骤：

（1）       gsl_monte_plain_state* gsl_monte_plain_alloc(size_t dim);

               申请一个计算dim重积分的gsl_monte_plain_state类型控制变量

（2）       int gsl_monte_plain_init(gsl_monte_plain_state* state);

               初始化控制变量

（3）       int gsl_monte_plain_integrate (const gsl_monte_function * f,

                           const double xl[], const double xu[],

                           const size_t dim,

                           const size_t calls,

                           gsl_rng * r,

                           gsl_monte_plain_state * state,

                           double *result, double *abserr);

        计算monte carlo积分，参数意义如下：

         ：f->被积分函数

         ：xl->积分下限数组，共有dim个

         ：xu->积分上限数组，共有dim个

         ：dim->积分重数

         ：calls->目标函数被调用次数，比如500000

         ：r->随机数产生器，比如r= gsl_rng_default

         ：state->状态空间，由gsl_monte_plain_alloc申请

         ：result->计算结果

         ：abserr->绝对误差

例子：计算多重积分

#include < iostream >
#include < time.h >
#include < gsl / gsl_rng.h >
#include < gsl / gsl_integration.h >
#include < gsl / gsl_monte.h >
#include < gsl / gsl_monte_miser.h >
#include < gsl / gsl_monte_plain.h >
#include < gsl / gsl_monte_vegas.h >
using namespace std;
 // visualsan@yahoo.cn SAN NUAA 2011.3.23
struct my_par
{
 double a,b,c;
};
 double my_f( double * x_array, size_t dim, void * params )
{
my_par * mp = (my_par * ) params ;
 if (dim != 2 )
{
fprintf(stderr, " dim!=2 " );
abort();
}
 double v = x_array[ 0 ] * x_array[ 0 ] * mp -> a +
x_array[ 0 ] * x_array[ 1 ] * mp -> b +
x_array[ 1 ] * x_array[ 1 ] * mp -> c;

 return exp( - 1 / v) / v;
}

 void test_monte()
{
my_par par = { 3 , 2 , 1 };
gsl_monte_function f;
f.dim = 2 ;
f.f = my_f;
f. params =& par;

 int calls = 500000 ;
 double xl[] = { 0 , 0 },xu[] = { 3.14 , 3.14 };
 double result,er;
gsl_monte_plain_state * ps = gsl_monte_plain_alloc( 2 );
 const gsl_rng_type * tp = gsl_rng_minstd;
gsl_rng * pr = gsl_rng_alloc(tp);
gsl_monte_plain_init(ps);
gsl_monte_plain_integrate( & f,
xl,xu, 2 ,calls,
pr,ps, & result, & er);


cout << result << endl; // ->0.913118
cout << er << endl; // ->0.0011711
/* matlab cmd:
fun = @(x,y) exp(-1./(3.*x.*x+2.*x.*y+y.*y))./( 3.*x.*x+2.*x.*y+y.*y );
quad2d(fun,0,3.14,0,3.14)
ans=0.913889482268682
 */

}
 void main()
{
test_monte();
} 


2.    miser monte carlo

     包含头文件：#include <gsl/gsl_monte_ miser.h>

  1. gsl_monte_miser_state* gsl_monte_miser_alloc(size_t dim);申请计算dim重积分的控制变量
  2. int gsl_monte_miser_init(gsl_monte_miser_state* state);初始化积分控制变量
  3. 积分函数：

        int gsl_monte_miser_integrate(gsl_monte_function * f,

                                         const double xl[], const double xh[],

                                         size_t dim, size_t calls,

                                         gsl_rng *r,

                                         gsl_monte_miser_state* state,

                                         double *result, double *abserr);

    4.void gsl_monte_miser_free(gsl_monte_miser_state* state);释放控制变量

miser carlo算法有几个参数需要设置，这些参数定义在gsl_monte_miser_state。

estimate_frac：这个参数在方差计算的递归调用过程中，指明每次函数调用次数占可用函数次数的百分比，默认值为0.1

min_calls：这个参数指明每次方差评估过程中函数调用的最小次数，如果申请函数调用次数n小于min_calls，则n=min_calls。min_calls默认值为16*dim.

min_calls_per_bisection：对分步骤中，函数调用的最小次数，默认值16* min_calls

alpha：指定两个子区域的方差如何组合，默认值为2

Dither：引入变异率，打破被积函数在积分区间的对称性。默认为0，当需要引入时，建议取值0.1

3.vegasmonte carlo

包含头文件gsl_monte_vegas.h

使用步骤：

1.）gsl_monte_vegas_state* gsl_monte_vegas_alloc(size_t dim);申请一个dim重积分的控制变量

2.）int gsl_monte_vegas_init(gsl_monte_vegas_state* state);初始化控制变量

3.）int gsl_monte_vegas_integrate(gsl_monte_function * f,

                                         double xl[], double xu[],

                                         size_t dim, size_t calls,

                                         gsl_rng * r,

                                         gsl_monte_vegas_state *state,

                                         double* result, double* abserr);

4.）void gsl_monte_vegas_free (gsl_monte_vegas_state* state);释放控制变量

下面这个例子是GSL手册上的例子：

计算积分

#include < iostream >
#include < time.h >
#include < gsl / gsl_math.h >
#include < gsl / gsl_monte_miser.h >
#include < gsl / gsl_monte_plain.h >
#include < gsl / gsl_monte_vegas.h >


using namespace std;
 // e: visualsan@yahoo.cn tr:SAN ad:NUAA t:2011
const double exact = 1.3932039296 ;
 double my_f( double * k, size_t dim, void * params )
{
 double A = 1.0 / ( M_PI * M_PI * M_PI );
 return A / ( 1.0 - cos(k[ 0 ]) * cos(k[ 1 ]) * cos(k[ 2 ]) );
}
 void disp_r( char * t, double r, double er, double ct)
{
printf( " %s----------------------t=%.1fms \n " ,t,ct);
printf( " result=%.6f\n " ,r);
printf( " sigma=%.6f\n " ,er);
printf( " exact=%.6f\n " ,exact);
printf( " error=%.6f = %.1g sigma\n " ,r - exact,fabs(exact - r) / er);
}
 void main()
{
cout << " monte carlo sample (visualsan@yahoo.cn tr:SAN )\n " ;
clock_t t1;
gsl_monte_function f;
f.dim = 3 ;
f.f = my_f;


 int calls = 500000 ;
 double xl[] = { 0 , 0 , 0 },xu[] = {M_PI,M_PI,M_PI};
 double result,er;

gsl_rng_env_setup();
 const gsl_rng_type * tp = gsl_rng_minstd;
gsl_rng * pr = gsl_rng_alloc(tp);



 // plain
t1 = clock();
gsl_monte_plain_state * ps = gsl_monte_plain_alloc( 3 );
gsl_monte_plain_init(ps);
gsl_monte_plain_integrate( & f,
xl,xu, 3 ,calls,
pr,ps, & result, & er);
disp_r( " plain " ,result,er, ( double )(clock() - t1) );
gsl_monte_plain_free(ps);

 // miser
t1 = clock();
gsl_monte_miser_state * pm = gsl_monte_miser_alloc( 3 );
gsl_monte_miser_init(pm);
gsl_monte_miser_integrate( & f,
xl,xu, 3 ,calls,
pr,pm, & result, & er);
disp_r( " miser " ,result,er,( double )(clock() - t1));
gsl_monte_miser_free(pm);

 // vegas
t1 = clock();
gsl_monte_vegas_state * pv = gsl_monte_vegas_alloc( 3 );
gsl_monte_vegas_init(pv);
gsl_monte_vegas_integrate( & f,
xl,xu, 3 , 10000 ,
pr,pv, & result, & er);
disp_r( " vegas warm-up " ,result,er,( double )(clock() - t1));

printf( " convering......\n " );
t1 = clock();
 do
{
gsl_monte_vegas_integrate( & f,
xl,xu, 3 ,calls / 5 ,
pr,pv, & result, & er);
printf( " result=%.6f,sigma=%.6f,chisq/dof=%.1f\n " ,result,er,pv -> chisq);
} while (fabs(pv -> chisq - 1.0 ) > 0.5 );
gsl_monte_vegas_free(pv);
disp_r( " vegas final " ,result,er,( double )(clock() - t1));

}

-------------------------------------------------------------------------------------------------

                                                                      SAN VISUALSAN@YAHOO.CN  2011.3.23  NUAA

转载于:https://www.cnblogs.com/JustHaveFun-SAN/archive/2011/03/23/san.html

weixin_30746117
关注
  • 0
    点赞
  • 踩
  • 0
    收藏
    觉得还不错? 一键收藏
  • 0
    评论
Cuba.jl：具有四个独立算法的多维数值积分库：Vegas，Suave，Divonne和Cuhre
02-04
Cuba.jl：具有四个独立算法的多维数值积分库：Vegas，Suave，Divonne和Cuhre
GSL----积分部分（翻译）
weixin_30394333的博客
03-21 368
GSL----积分部分（翻译） TR:SAN EMIAL:VISUALSAN@YAHOO.CN NUAA 2011.3.22 --------------------------------------------------------------- 每个算法都是计算积分表达式的近似值： 其中w(x)是权函数，一般积分时取w(x)=1.通过提供绝对精度epsabs和相对精度ep...
蒙特卡洛积分
Knight_Ren的博客
11-28 2119
描述 蒙特卡洛积分是利用随机试验求函数的数值积分的方法。请考虑定义一个蒙特卡洛函数mt_int，通过生成平面区域[x1,x2]*[y1,y2]里的随机点数据，检查其位于被积函数f的值之下的比率，得到函数积分的近似值。 本题目中被积函数为f(x)=(x2+x+ex)/(x3-x+1) 用蒙特卡洛积分求该函数与x轴的面积。 输入 N+1行 第一行输入整数N，表示接下来会有N行输入 接下来每行输入一个浮点数x。( 0 <= x <= 7.0) 输出 被积函数f在区间[0,x]内曲线与
蒙特卡洛（Monte Carlo）积分的入门
CS笔记本
05-12 2124
转自 http://www.opengpu.org/forum.php?mod=viewthread&tid=198 概率密度描述了一个随机变量的值在区域的概率分布。 概率密度函数（probability density function）为描述随机变量的概率分布。亦为累积分布函数（CDF）之导函数 定义 对于一维实随机变量X，任何一个满足下列条件的函数fX(x)都
概率论和蒙特卡洛积分
averagePerson的博客
07-19 913
离散随机变量连续随机变量与PDF。
GSL中的数值积分
yxmlhc的专栏
09-28 996
数值积分 本章描述了对一维函数进行数值积分(求积)的程序。常用函数有自适应和非自适应集成的例程，具体情况有专门的例程。这些包括在无限和半无限范围内的积分，奇异积分，包括对数奇点，柯西主值的计算和振荡积分。本库重新实现了在QUADPACK中使用的算法，数值积分包由Piessens、de donker - kapenga、Ueberhuber和Kahaner编写。在Netlib上可以找到Fortran版的QUADPACK代码。包括具有高精度系数的非自适应、定阶Gauss-Legendre积分例程，以及IQPA
gsl-ref.rar_gsl_gsl 函数_gsl-1.8
09-20
《gsl-ref.rar_gsl_gsl 函数_gsl-1.8》是一个关于GNU Scientific Library (GSL)的参考资料，特别关注的是GSL函数的详细描述，适用于版本1.8。GSL是一个开源的数学库，它提供了广泛的数值算法，包括线性代数、微积分...
GSL-2.4.tar.gz
08-11
首先，你需要从官方仓库或者镜像站点下载 GSL-2.4 的源代码包，文件名为 "GSL-2.4.tar.gz"。使用 `wget` 或者浏览器下载后，使用 `tar` 命令解压： ``` wget http://ftp.gnu.org/gnu/gsl/GSL-2.4.tar.gz tar -...
gsl-1.15.tar.gz
02-28
标题 "gsl-1.15.tar.gz" 指的是GNU Scientific Library（GSL）的版本1.15的源代码压缩包。GSL是一个广泛使用的开源库，旨在为C语言提供各种科学和工程计算所需的算法。这个库涵盖了线性代数、数值积分、特殊函数、...
gsl-1.16 vc2012
02-28
标题 "gsl-1.16 vc2012" 指的是GNU Scientific Library (GSL)的1.16版本，针对Microsoft Visual Studio 2012编译环境进行了编译。GSL是一个开源的C语言库，提供广泛的数学、统计和数值计算功能，适用于科学和工程...
gsl-1.5-gsl-1.6.patch.gz_正态分布
07-14
可以生成十几种各类分布的随机数，如泊松分布，正态分布，指数分布等等
GSL积分运算
Augusdi的专栏
08-15 2789
#include #include #include #pragma comment(lib, "libgsl_d.lib") #pragma comment(lib, "libgslcblas_d.lib") double f (double x, void * params) { double alpha = *(double *) params; double f = lo
GSL----积分部分（翻译+精简）
Augmented Reality & Computer Graphics
03-21 3411
GSL----积分部分（翻译+精简） SAN每个算法都是计算积分表达式的近似值：其中w(x)是权函数，一般积分时取w(x)=1.通过提供绝对精度epsabs和相对精度epsrel，达到预期的计算精度满足：注意这是不是一个同时满足的约束条件，而是当满足绝对误差或者相对误差或者都满足时计算结束。可以令epsabs=0或epsrel=0,这表明你想计算积分的精确值，此时gsl可能会因为条件太苛刻而计算失败，但是通常情况下gsl会给出一个尽可能满意的结果。积分函数命名中满足下列简称规则：Q:  quadrature
GSL数值积分例子
Augusdi的专栏
08-14 3307
//数值计算实验 数值积分 #include #include #include #include using namespace std; //被积函数 double f(double x) { //为便于调试，先弄个有原函数的 y = x^2 + x^3 - 2*x^4 return x*x + x*x*x - 2*x*x*x*x; } //被积函数 给gsl用的
建模：蒙特卡洛法求多重积分
生姜用户的博客
08-24 1846
蒙特卡洛法求多重积分-Cpp 0~1随机分布的点列 C++类的写法： class RND { private: int R; public: RND(int r = 1) { R = r; } double rnd1() { int m; double s, u, v, p; s = 65536.0; u = 2053.0; v = 13849.0; ...
GSL科学计算库——计算高斯-勒让德积分
ouening的博客
12-28 3664
假设计算下列积分： ∫0πexcos(x)dx\int_0^\pi e^xcos(x)dx∫0π​excos(x)dx 根据高斯勒让德积分，先将积分区间变换为[−1,1][-1,1][−1,1]，使用线性变换： x=π−02t+π+02x=\frac{\pi-0}{2}t+\frac{\pi+0}{2}x=2π−0​t+2π+0​，将原积分变为 ∫0πexcos(x)dx=−π2eπ2t∫−11e...
GSL蒙特卡洛积分
Augusdi的专栏
08-15 2120
#define GSL_DLL #include #include #include #include #include #include #pragma comment(lib, "libgsl_dll_d.lib") #pragma comment(lib, "libgslcblas_dll_d.lib") /* Computation of the integral,
C/C++ 之 GSL 数学运算库使用笔记
最新发布
一些笔记
03-31 2076
本文介绍了GSL 数学运算库的安装、编译与使用教程，小白友好，放心食用~
蒙特卡洛（Monte-Carlo）计算定积分【Python】
热门推荐
Hren0412的博客
10-11 1万+
1、蒙特卡洛计算定积分【投点法】： 以求解第一小题（一重积分）为例，如下图所示，蓝色曲线为被积函数f(x)在[0, π/2]区间上的图像，若要求其定积分，其实就是求曲线下方的面积。这时我们可以用一个比较容易算得面积的矩型罩在函数的积分区间上（假设其面积为S）。然后随机地向这个矩形框里面投点，统计出蓝色曲线下方的点的数量占所有点（红框内）数量的比例为R，那么就可以据此估算出被积函数从0到π/2的定...
Ubuntu环境下详细安装ns-3教程
sudo apt-get install gsl-bin libgsl0-dev libgsl0ldbl ``` 为了进行词法和语法分析，安装Flex和Bison： ``` sudo apt-get install flex bison ``` 确保C++编译器和相关工具也已安装： ``` sudo apt-get install ...
评论
添加红包

请填写红包祝福语或标题

个

红包个数最小为10个

元

红包金额最低5元

当前余额3.43元 前往充值 >
需支付：10.00元
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付元
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明：

1.余额是钱包充值的虚拟货币，按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载，可以购买VIP、付费专栏及课程。

余额充值