#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; }