Вычисление интегралов методом Монте-Карло

  bord: tBorders; {массив границ интегрирования}

 

begin

  randomize;

  writeln('Выберите распределение случайной величины:');

  writeln('1 - Равномерное распределение');

  writeln('2 - Нормальное распределение');

  writeln('3 - Экспоненциальное распределение');

  writeln('4 - Бета распределение');

  writeln('5 - Гамма-распределение');

  readln(r);

  Readln(mul);

  for i := mul downto 1 do begin

    write('нижняя  граница: ', mul-i+1, ' = ');

    Readln(bord[i].bottom); {i-я нижняя граница интегрирования}

    write('верхняя граница ', mul-i+1, ' = ');

    Readln(bord[i].top); {i-я верхняя граница интегрирования}

  end;

  write('Количество испытаний = ');

  readln(n);

  case r of

    1:  {вычисление интеграла по равномерному  распределению}

      begin

        sum := 0;

        sum2 := 0;

        for i := 1 to n do

        begin

          pMult := 1;

          for j := 1 to mul do begin

            ksi[j]:= bord[j].bottom + (bord[j].top -               bord[j].bottom) * random;

            pMult := pMult * pMultipleUniformDist(j, bord);

          end;

          val := f(ksi, bord, mul) / pMult ;

          sum := sum + val;

          sum2 := sum2 + val * val;

        end;

        writeln('Результат: ', sum / n: 0: 8);

        disp := (sum2 - sum * sum / n) / n;

        writeln('Ошибка: ', 3 * sqrt(disp / n): 0: 8);

        readln;

      end;

 

    2: {вычисление интеграла по нормальному распределению}

      begin

        sum := 0;

        sum2 := 0;

        for i := 1 to n do

        begin

          pMult := 1;

          for j := 1 to mul do begin

            ksi[j] := randomValueNormDist;

            pMult := pMult * pNormDist(ksi[j]);

          end;

          val := f(ksi, bord, mul) / pMult;

          sum := sum + val;

          sum2 := sum2 + val * val;

        end;

        writeln('Результат: ', sum / n: 0: 8);

        disp := (sum2 - sum * sum / n) / n;

        writeln('Ошибка: ', 3 * sqrt(disp / n): 0: 8);

        readln;

      end;

 

    3:  {вычисление интеграла по экспоненциальному распределению}

      begin

       write('L = ');

        readln(L);

        sum := 0;

        sum2 := 0;

        for i := 1 to n do

        begin

          pMult := 1;

          for j := 1 to mul do begin

            alpha := random;

            ksi[j] := randomValueExpDist(L);

            pMult := pMult * pExpDist(ksi[j], L);

          end;

          val := f(ksi, bord, mul) / pMult;

          sum := sum + val;

          sum2 := sum2 + val * val;

        end;

        writeln('Результат: ', sum / n: 0: 8);

        disp := (sum2 - sum * sum / n) / n;

        writeln('Ошибка: ', 3 * sqrt(disp / n): 0: 8);

        readln;

      end;

 

    4:  {вычисление интеграла по бета-распределению}

      begin

        write('p = ');

        readln(p);

        write('m = ');

        readln(m);

        sum := 0;

        sum2 := 0;

        for i := 1 to n do

        begin

          pMult := 1;

          for j := 1 to mul do begin

            alpha := random;

            ksi[j] := randomValueBetaDist(m, p);

            pMult := pMult * pBetaDist(ksi[j], p, m);

          end;

          val := f(ksi, bord, mul) / pMult;

          sum := sum + val;

          sum2 := sum2 + val * val;

        end;

        writeln('Результат: ', sum / n: 0: 8);

        disp := (sum2 - sum * sum / n) / n;

        writeln('Ошибка: ', 3 * sqrt(disp / n): 0: 8);

        readln;

      end;

 

    5:  {вычисление интеграла по гамма-распределению}

      begin

        write('W= ');

        readln(W);

        write('B= ');

        readln(B);

        sum := 0;

        sum2 := 0;

        for i := 1 to n do

        begin

          pMult := 1;

          for j := 1 to mul do begin

            ksi[j] := randomValueGammaDist(W, B);

            pMult := pMult * pGammaDist(ksi[j], W, B);

          end;

          val := f(ksi, bord, mul) / pMult;

         sum := sum + val;

          sum2 := sum2 + val * val;

       end;

        writeln('Результат: ', sum / n: 0: 8);

        disp := (sum2 - sum * sum / n) / n;

        writeln('Ошибка: ', 3 * sqrt(disp / n): 0: 8);

        readln;

      end;

   end;

end;

 

Пример работы программы:

 

Вычислим интеграл

 

Точное значение интеграла:  0,25

 

 

Результаты расчётов программы:

 

	Приближённое значение интеграла
	Равномерное распределение		Нормальное распределение		Экспоненциальное распределение		Гамма-распределение		Бета-Распределение
Количество испытаний	Результат	Ошибка	Результат	Ошибка	Результат	Ошибка	Результат	Ошибка	Результат	Ошибка
n = 10000	0,24817	0,00661	0,24441	0,03325	0,24338	0,01992	0,24791	0,00655	0,25613	0,02060
n = 100000	0,25088	0,00209	0,25316	0,01058	0,25070	0,00638	0,25013	0,00208	0,25177	0,00641
n = 1000000	0,25029	0,00066	0,25021	0,00332	0,25041	0,00202	0,25025	0,00066	0,25126	0,00202

 

 

 

Приложение 2

Программа для решения интегрального уравнения Вольтерра методом Монте-Карло (Delphi)

 

 

var

  n: longint; {объём выборки}

  t: double; {верхняя граница интегрирования}

  sum, sum2: double; {суммы вычисленных значений и  квадратов этих значений}

  disp: double; {дисперсия}

  t0, t1, q: double;

begin

  write('top border(t) = ');

  readln(t);

  sum := 0;

  sum2 := 0;

  for i := 1 to n do begin

    sum := sum + y(t);

    sum2 := sum2 + y(t)*y(t);

    t0 := t;

    t1 := Random;

    q := 1;

    while (t1 < t0) do begin

      q := q * a(t1);

      val := y(t1)*q;

    sum := sum + val;

      sum2 := sum2 + val * val;

      t0 := t1;

      t1 := Random;

    end;

  end;

  writeln('Результат: ', sum/n : 0: 8);

  disp := (sum2 - sum * sum / n) / n;

  writeln('Ошибка: ', 3 * sqrt(disp / n): 0: 8);

  readln;

 end;

 

Пример работы программы:

 

Вычислим интегральное уравнение  Вольтерра  при

 

Точное значение уравнения: 

 

 

Объём выборки	Приближённое решение  уравнения	Ошибка
n = 100000	0,35856	0,00227
n = 1000000	0,35995	0,00082
n = 10000000	0,36006	0,00008

Результаты расчётов программы:

 

 

 

Заключение

 

В процессе исследований было выяснено, что метод Монте-Карло имеет несколько важных преимуществ:

1. Метод прост в реализации;
2. Он легко применим без предварительного анализа решаемой задачи;
3. Неплохо справляется с вычислением интегралов большой кратности, которые не вычислить другими численными методами.

Также были замечены очевидные недостатки метода:

1. По результатам работы программ можно заметить, что для получения хорошего результата нужно проводить сотни тысяч испытаний;
2. Без наличия датчика псевдослучайных чисел метод неприменим;
3. При увеличении количества испытаний погрешность расчётов убывает медленно.

 

 

Литература

1. С. М. Ермаков, С. А. Михайлов «Статистическое моделирование», 1982г.
2. С. М. Ермаков «Метод Монте-Карло в вычислительной математике», 2009г.