Функция Monte_Carlo делает вещи более сложными, чем они должны быть.Для интеграции одномерной функции все, что нам нужно сделать, это сэмплировать значение функции несколько раз в пределах области, в которую мы интегрируем:
#include <random>
double Monte_Carlo(double Func(double), double xmin, double xmax, int N)
{
// This is the distribution we're using to generate inputs
auto x_dist = std::uniform_real_distribution<>(xmin, xmax);
// This is the random number generator itself
auto rng = std::default_random_engine();
// Calculate the total of N random samples
double total = 0.0;
for(int i = 0; i < N; i++) {
double x = x_dist(rng); // Generate a value
total += Func(x);
}
// Return the size of the interval times the total,
// divided by the number of samples
return (xmax - xmin) * total / N;
}
Если мы запустим этот код с N = 1000
, мы получаем целое значение 9.20569
, которое довольно близко к точному ответу (9.33333...
).
// It's much more efficent to use x*x instead of pow
double func(double x) { return x * x + 1; }
int main()
{
cout << "Integral value is: " << Monte_Carlo(func, -2, 2, 1000) << endl;
getchar(); // Pause until the user presses enter
return 0;
}
Мы также можем попробовать несколько значений N
, чтобы программа показывала, как она сходится.Следующая программа вычисляет интеграл с N, являющимися степенями 2
от 0
до 30
#include <iostream>
#include <cmath>
#include <random>
using namespace std;
double func(double x) { return x*x + 1; }
double Monte_Carlo(double Func(double), double xmin, double xmax, int N) {
auto x_dist = std::uniform_real_distribution<>(xmin, xmax);
auto rng = std::default_random_engine();
double total = 0.0;
for(int i = 0; i < N; i++) {
double x = x_dist(rng); // Generate a value
total += Func(x);
}
return (xmax - xmin) * total / N;
}
int main() {
int N = 1;
for(int i = 0; i < 31; i++) {
std::cout << "N = " << N << "\t\tintegral = " << Monte_Carlo(func, -2, 2, N) << endl;
N *= 2; // Double N
}
}
Выходные данные показывают, что метод Монте-Карло действительно сходится:
N = 1 integral = 12.6889
N = 2 integral = 8.39917
N = 4 integral = 7.97521
N = 8 integral = 9.24233
N = 16 integral = 9.75632
N = 32 integral = 9.87064
N = 64 integral = 9.46945
N = 128 integral = 9.27281
N = 256 integral = 9.27395
N = 512 integral = 9.17546
N = 1024 integral = 9.19097
N = 2048 integral = 9.26203
N = 4096 integral = 9.37979
N = 8192 integral = 9.36167
N = 16384 integral = 9.28918
N = 32768 integral = 9.29766
N = 65536 integral = 9.31101
N = 131072 integral = 9.3227
N = 262144 integral = 9.32588
N = 524288 integral = 9.32805
N = 1048576 integral = 9.32726
N = 2097152 integral = 9.32722
N = 4194304 integral = 9.331
N = 8388608 integral = 9.33082
N = 16777216 integral = 9.33174
N = 33554432 integral = 9.33164
N = 67108864 integral = 9.33303
N = 134217728 integral = 9.33283
N = 268435456 integral = 9.33327
N = 536870912 integral = 9.33325
N = 1073741824 integral = 9.33333