С++ решение интеграла методом Монте-Карло (нужно проверить что не так в коде)

1. Зайди на сайт http://yotx.ru и набери свою функцию: её максимум совсем не на конце отрезка, а в районе 0.525. Потому с и d ты считаешь неправильно.

2. a и b - вещественные переменные.

Проще всего: double a = 0.0, b = 1.0, с = 0.0, d = 0.25; - безо всяких вычислений.

3. Зачем тебе массивы x и y? Тебе нужны только текущие значения x и y.

4. Про rand тебе сказали - зря не послушал. Правильное вычисление:
x = a + (b - a) * rand() / RAND_MAX;
y = c + (d - c) * rand() / RAND_MAX;
Но т. к. a, b, с и d известны, то формулы упрощаются до:
x = rand() * 1.0 / RAND_MAX; // ОБЯЗАТЕЛЬНО вещественное деление
y = rand() * 0.25 / RAND_MAX;

5. Не y[i] < func(x[i], а y <= func(x): равно - это тоже в границах интеграла.

В результате всё вычисление сокращается до:

long k = 0, q;
cout << "q = ";
cin >> q;
srand(time(NULL));
for (long i = 0; i < q; ++i) {
if (rand() * 0.25 / RAND_MAX <= func(rand() * 1.0 / RAND_MAX)) { ++k; }
}
cout << "I = " << k * 0.25 / q;

Я вообще по другому сделал. у не может превышать 0,5, поэтому х задаётся в диапазоне [0;1], а у в диапазоне [0;0.5]. А дальше - обычный алгоритм:

double r()

{ return rand()%1000001*1e-7; }

int main()

{ double x,y,s,S=0.185376068; int k,l,n; srand(time(NULL)); for(;;) { cout <<

"n = ?\b"; cin >> n; l=0; for(k=0; k < n; k++) { x=r(); y=0.5*r(); if (2*y*(1+x*(1+x)) <= sin(2*x)) l++; } s=0.5*l/n; cout << s << ' ' << s/S-1 << '\n'; } }

Зачем считать синус и косинус, когда можно считать только один раз синус, что сразу даёт выигрыш в оперативности раза в два? А результаты такие (несколько раз количество испытаний n я задавал по миллиону):

0.1857535 (относительная ошибка 2.03е-3)

0.185589 (ошибка 1.149е-3)

0.185386 (ошибка 5.897е-5, неожиданный рекорд точности!)

В общем и целом ошибка получается порядка ~¹/√n и плохо сглаживается увеличением количества испытаний, поэтому увеличивать n, например, до миллиарда - бессмысленно!..

rand генерирует целые числа, а у тебя диапазон от 0 до 1