#include <stdio.h> //printf
#include <time.h> //用时间作为随机种子
#include <stdlib.h> //rand和srand
#include <math.h> //数学函数
//宏定义
#define max(a, b) ((a) > (b) ? (a) : (b))
#define min(a, b) ((a) > (b) ? (b) : (a))
inline double random()
{
return rand() % 32768 / 32767.0;
}
//Monte carro求给定函数的定积分
//极限:[-1, 1.414213562373095]
inline double f1(double x, double xmin, double xmax)
{
if (x >= xmin && x <= xmax)
return x * x * sqrt(1 + x * x * x) - 1;
}
//极限:[0, 0.693147180559945]
inline double f2(double x, double xmin, double xmax)
{
if (x >= xmin && x <= xmax)
return x * log(1 + x);
}
//极限:[0, 1]
inline double f3(double x, double xmin, double xmax)
{
if (x >= xmin && x <= xmax)
return sqrt(1 - x * x);
}
//极限:[-0.35969769413186028259906339255702, 0.1]
inline double f4(double x, double xmin, double xmax)
{
if (x >= xmin && x <= xmax)
return cos(x) - 0.9;
}
//定义函数指针类型
typedef double (*Func)(double x, double xmin, double xmax);
//定义结构体存储函数的上下限
typedef struct Limit {
double lVal, hVal;
//C++结构体构造函数
Limit(double lv, double rv)
{
lVal = lv;
hVal = rv;
}
} Limit;
//取得函数的上下极限——其精确度影响很大---
//由于用于积分,故若y>0,底限0;若y<0,则高限0;若y穿越x轴,则不修改默认极限
Limit getLimit(Func func, double xmin, double xmax)
{
//边界极限
double t1 = func(xmin, xmin, xmax), t2 = func(xmax, xmin, xmax);
Limit limit(min(t1, t2), max(t1, t2));
double delta = 1e-3, lastX, lastY, y;
int i = 0;
for(double x = xmin; x <= xmax; x += delta, i++)
{
y = func(x, xmin, xmax);
if (y < limit.lVal)
limit.lVal = y;
if (y > limit.hVal)
limit.hVal = y;
if (i == 0)
{
lastX = x;
lastY = y;
}
else
{
double tx = (lastX + x) / 2;
double ty = func(tx, xmin, xmax);
if (fabs(lastY) > fabs(y))
{
if (fabs(lastY) < fabs(ty))
{
lastX = tx;
lastY = ty;
if (ty < 0 && limit.lVal > ty)
limit.lVal = ty;
else if (ty >= 0 && limit.hVal <= ty)
limit.hVal = ty;
}
}
if (fabs(lastY) < fabs(y))
{
if (fabs(y) < fabs(ty))
{
lastX = tx;
lastY = ty;
if (ty < 0 && limit.lVal > ty)
limit.lVal = ty;
else if (ty >= 0 && limit.hVal <= ty)
limit.hVal = ty;
}
}
}
}
if (limit.lVal > 0)
limit.lVal = 0;
if (limit.hVal < 0)
limit.hVal = 0;
return limit;
}
//积分函数——func:函数指针
double cumulate(Func func, int N, double xmin, double xmax)
{
Limit limit = getLimit(func, xmin, xmax);
double x, y, y0;
int cntA = 0, cntB = 0;
for(int i = 0; i < N; i++)
{
//随机在平面上投掷筛子,位置(x, y)
x = xmin + (xmax - xmin) * random();
y = limit.lVal + (limit.hVal - limit.lVal) * random();
//调用函数取x点位置的标准y值
y0 = func(x, xmin, xmax);
//若y和yr都在x轴同侧
if ( y0 >= 0 && y >= 0 && y <= y0) //正积分计数
cntA++;
else if (y0 <= 0 && y <= 0 && y >= y0) //负积分计数
cntB++;
}
printf("Limit(%.15f, %.15f\n", limit.lVal, limit.hVal);
return (cntA - cntB) * (xmax - xmin) * (limit.hVal - limit.lVal) / N;
}
int main()
{
double rst;
//每个函数随机投筛子N次
int N = 1e7;
//函数指针数组
Func fs[] = {f1, f2, f3, f4};
//初始化随机种子
srand (time(NULL));
//经Matlab对上述函数积分得到的结果
double ans[] = {-0.593682861167513, 0.25, 3.141592653589793/4, -0.058529015192103493347497678369701};
//分别调用4个函数
for(int i = 0; i < 4; i++)
{
rst = cumulate(fs[i], N, 0, 1);
printf("标准:%.15f, 程序结果:%.15f\n", ans[i], rst);
}
return 0;
}
蒙特卡洛(Monte Carlo)自动求函数积分的C++算法实现与测试
最新推荐文章于 2022-01-12 22:05:04 发布