Метод Монте-Карло - какие-то жуткие значения

[byulent]

Пытаюсь написать программу для интегрирования методом Монте-Карло.
По задаче:

$$\int\limits^1_0x^{-1/4}(1-x)^{-3/4}dx=\frac{\pi}{\sin({3\pi/4})}$$

Мало того, что WolframAlpha выдаёт, что точное значение в 2 раза больше, так ещё и в самих вычислениях какая-то хрень творится.
Пусть $\xi_i\in[0;1]$ - это у нас

Используется синтаксис C++

(double)(rand())/RAND_MAX;

,
а $\rho(x)=1$. Тогда $x_i=\xi_i+a$.
Написал такой код:

Используется синтаксис C++

for (int i=0; i<N; i++){
        x[i]=r()+a;
        //cout<<x[i]<<" "<<f(x[i])<<endl;
        sum+=f(x[i]);
    }
    //cout<<sum<<endl;
    I1=sum/N;
    for (int i=0; i<N; i++) sq+=pow(f(x[i])-I1,2);
    disp1=sqrt(sq/(N-1));
    cout<<"При rho=1:"<<endl<<"I^="<<I1<<", sigma^="<<disp1<<endl;

Результаты:
$N=10: \hat I=2,915; \hat\sigma=1,68875$
$N=100: \hat I=3,32821; \hat\sigma=3,17415...$
Понял, что фигня какая-то, попробовал сменить $\rho: \rho(x)=2(x+2)/3$.

Используется синтаксис C++

    for (int i=0; i<N; i++){
        x[i]=-2+sqrt(4+2*r()/C);
        cout<<x[i]<<" "<<f(x[i])/(C*(x[i]+2))<<endl;
        sum+=f(x[i])/(C*(x[i]+2));
    }
    I2=sum/N;
    for (int i=0; i<N; i++) sq+=pow(f(x[i])/(C*(x[i]+2))-I2,2);
    disp2=sqrt(sq/(N-1));
    cout<<"При rho=x+2:"<<endl<<"I^="<<I2<<", sigma^="<<disp2<<endl;

Лучше не стало:

$N=10: \hat I=3,60127; \hat\sigma=10,1102$
Почему у интеграла получаются такие значения и какая-то жуткая дисперсия? Что здесь не так?

[Aritaborian]

byulent в сообщении #1129593 писал(а):

Мало того, что WolframAlpha выдаёт, что точное значение в 2 раза больше

Всё она верно выдаёт.

[svv]

Сходимость очень плохая из-за того, что область между графиком $f(x)=x^{-1/4}(1-x)^{-3/4}$ и осью абсцисс неограниченная (потому что сама функция неограничена). Даже при миллионе точек точность неприемлемая. Чтобы улучшить сходимость, можно взять такую функцию $f_0(x)$, чтобы интеграл $\int\limits_0^1 f_0(x)\;dx$ легко брался, а разность $f(x)-f_0(x)$ была уже ограниченной на $[0,1]$. Например,
$f_0(x)=x^{-1/4}+(1-x)^{-3/4}$
Ну, и искать методом Монте-Карло интеграл
$\int\limits_0^1 (f(x)-f_0(x))\;dx$,
прибавив затем найденный аналитически $\int\limits_0^1 f_0(x)\;dx$.
Я проверил, это помогает, хотя точек всё равно надо брать много.

[byulent]

Aritaborian в сообщении #1129597 писал(а):

Всё она верно выдаёт.

Да, извиняюсь, забыл, что $2/\sqrt{2}=\sqrt{2}$.

-- 07.06.2016, 07:21 --

В коде ошибок точно нет?

Используется синтаксис C++

double f (double x){return pow(x,-0.25)*pow(1-x,-0.75);} //функция
double r (){return (double)(rand())/RAND_MAX;} //случайное число, равномерно распределённое на единичном отрезке

int main()
{
    /**
    * Метод Монте-Карло
    **/
    setlocale(LC_ALL,"rus");
    double I1, disp1, I2, disp2, I3, disp3, a, b, *x, sum=0, sq=0;
    int N;
    cout<<"Введите пределы интегрирования"<<endl;
    cin>>a>>b;
    cout<<"Введите количество испытаний"<<endl;
    cin>>N;
    x = new double[N];
    srand(time(0));
    //p=1
    for (int i=0; i<N; i++){
        x[i]=r()+a;
        //cout<<x[i]<<" "<<f(x[i])<<endl;
        sum+=f(x[i]);
    }
    //cout<<sum<<endl;
    I1=sum/N;
    for (int i=0; i<N; i++) sq+=pow(f(x[i])-I1,2);
    disp1=sqrt(sq/(N-1));
    cout<<"При rho=1:"<<endl<<"I^="<<I1<<", sigma^="<<disp1<<endl;
    //p=x+2
    sum=0;
    double C=1/(pow(b,2)/2+2*b-pow(a,2)/2-2*a);
    //cout<<C<<endl;
    for (int i=0; i<N; i++){
        x[i]=-2+sqrt(4+2*r()/C);
        //cout<<x[i]<<" "<<f(x[i])/(C*(x[i]+2))<<endl;
        sum+=f(x[i])/(C*(x[i]+2));
    }
    I2=sum/N;
    for (int i=0; i<N; i++) sq+=pow(f(x[i])/(C*(x[i]+2))-I2,2);
    disp2=sqrt(sq/(N-1));
    cout<<"При rho=x+2:"<<endl<<"I^="<<I2<<", sigma^="<<disp2<<endl;
    return 0;
}

[ИСН]

Код проверьте на какой-нибудь нормальной функции, без бесконечностей. Лучше на нескольких.