Паскаль. Основы программирования


           

с помощью генератора случайных чисел


          begin

            x := random*(b - a) + a;

            f := f + fx(x)

          end;

        I := (b - a)*f/n

      end;

Итак, с помощью генератора случайных чисел (random) вырабатывается случайное число из промежутка [a, b], находится значение функции в этой точке и суммируются ее абсолютные величины (учитывая, что функция может принимать и отрицательные значения) в переменную f. Приближенное значение интеграла вычисляется с помощью оператора I := (b - a)*f/n.

Теперь стоит очень болезненный вопрос оценки точности значения интеграла. Не вдаваясь в долгие математические рассуждения (о чем смотрите [19]), покажем неравенства, с помощью которых можно оценить точность вычисления интеграла.

Пусть I(f) - точное значение интеграла, а Sn(f) - его приближенное значение, тогда с вероятностью 0,997 и 0,99999 выполняются следующие соотношения

 и 


Здесь D(f) - дисперсия непрерывной случайной величины, которая вычисляется по формуле:



где сумма есть математическое ожидание случайной величины.

Приведем схему последовательного вычисления интеграла с заданной точностью eps. Последовательно, при n = 1, ... получаются случайные точки x и вычисляющая величины tn

, Sn , dn , пользуясь рекуррентными  соотношениями









и величину  
  или 


Начальные условия рекурсии n = 1,
 
,


Если оказалось, что
 <= eps или
 <=eps, то вычисления прекращаются и полагают, что приближенное значение интеграла равно Sn с вероятностью 0,997 или 0,99999 и точностью eps.

Изложенные математические соображения можно реализовать следующей процедурой:

{ Процедура оценки точности и вычисления интеграла }

   Procedure Monte_Karlo(a, b, eps : real; var

s : real);

      var

        t, d, f, dd  : real;

        n               : longint;

      begin

        n := 1;

        t := I(n, a, b);

        s := t;

        d :=0;

          repeat

            n := n + 1;

            t := t + I(n, a, b);

            s := t/n;

            d := d + (n/(n - 1))*sqr(I(n, a, b) - s);

            dd := d/(n - 1)


Содержание  Назад  Вперед