Интегрирование методом Монте-Карло

Question

Пытаюсь найти площадь sin(x)/x В границах x1 = 0.001, x2 = 11, y1 = -0.2, y2 = 1; Ответ почему-то не сходится с правильным. Уже всё перепроверил. Подскажите в чём ошибка. Правильный ответ S = 1.57, количество точек попавшие n = 2400. Мой ответ S= 2.6, количество точек попавшее n = 4000

#define _USE_MATH_DEFINES
#include <iostream>
#include <ctime>
#include <cmath>

//бросает точку
double random(double min, double max) {

    double r = (double)rand() / RAND_MAX;

    return min + r*(max - min);
}
//вычисляет Y
double dop(double x) {
    //задаём формулу
    double ydop = (sin(x))/x;

    return ydop;
}

double border(double min, double max) {
    return 1;
}


int main()
{
    srand(time(NULL));
    //задаём границы
    double x1 = 0.001, x2 = 11, y1 = -0.2, y2 = 1;
    //число бросков
    double n = 20000;

    int W = 0;
    double S = 0;
   
    double Yrand, Xrand, Ydop;
    for (int i = 0; i < n; i++) {
        Xrand = random(x1, x2);
        Yrand = random(y1, y2);
        Ydop = dop(Xrand);
        if (Ydop > Yrand && Ydop > 0 && Yrand > 0) {
            W++;
        }
        else if (Ydop < Yrand && Ydop < 0 && Yrand < 0) {
            W++;
        }
    }

    S = abs(x1 - x2) * abs(y1 - y2);

    double ot = (static_cast<double>(W) / n) * S;
    std::cout << W << std::endl;
    std::cout << ot;
}

Answer 1

Я бы делал так:

for (int i = 0; i < n; i++) {
    Xrand = random(x1, x2);
    Yrand = random(y1, y2);
    Ydop = dop(Xrand);
    if (Ydop > Yrand) W++;
}

S = abs(x1 - x2) * abs(y1 - y2);

double ot = (static_cast<double>(W) / n) * S + y1*(x2-x1);

Поскольку при интегрировании sin(x)/x надо не забывать, что то, что ниже оси абсцисс, идет со знаком минус...

Но и тут результат будет неверный, потому что функция опускается ниже -0.2, так что я бы еще исправил y1 = -0.3. Тогда результат становится похож на правду.

Второй вариант чуть сложнее — учет знаков (но y1 = -0.3 все равно надо оставить):

double Yrand, Xrand, Ydop;
for (int i = 0; i < n; i++) {
    Xrand = random(x1, x2);
    Yrand = random(y1, y2);
    Ydop = dop(Xrand);
    if (Ydop >= 0 && Ydop >= Yrand && Yrand >= 0) W++;
    else if (Ydop < 0 && Ydop <= Yrand && Yrand <= 0) W--;
}

S = abs(x1 - x2) * abs(y1 - y2);

double ot = (static_cast<double>(W) / n) * S;

БЛОГ НА HUSL

Интегрирование методом Монте-Карло

Ответы (1 шт):