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

Алгоритм мне прислал Nick Maslov. Функция вычисляет интеграл вида:

Мы будем считать, что область D вписана в квадрат [a_x , b_x ]*[a_y , b_y ], тогда полагаем:

и используем для приближенного вычисления интеграла формулу

где (x_i , y_i ) пара независимых равномерно распределенных случайных величин.

На вход функции подается параметры ax, ay, bx, by, которые определяют квадрат [a_x , b_x ]*[a_y , b_y ] на плоскости, в который вписана область D. В блок-схеме предполагается, что определена функция F(x,y), интеграл от которой необходимо вычислить, а так же функция D(x,y), которая принимает значение истина в случае, когда точка (x,y) принадлежит области интегрирования D и ложь в противном случае.

Параметр n опеределяет кол-во рассматриваемых точек и соответственно влияет на точность вычисления интеграла, но не следует брать его слишком большим иначе придется изменить алгоритм, так как может произойти переполнение при суммировании. Алгоритм нетрудно обобщить для вычисления интегралов от многомерных функций.

Если нашли ошибку в алгоритме - сообщите!

Реализации алгоритма на различных языках:

Блоксхемы:

Реализация алгоритма на AlgoPascal:

unit MonteCarloIntegrationUnit;

declaration
    //функция, задающая область интегрирования
    function D(X:Real; Y:Real):Boolean;

    //интегрируемая функция
    function F(X:Real; Y:Real):Real;

interface MonteCarloIntegration;

implementation

(*
Интегрирование методом Монте-Карло на области D, включенной
в прямоугольник [ax, bx]*[ay, by], с использованием n опорных
точек.

*)
function MonteCarloIntegration(
    const ax : Real;
    const bx : Real;
    const ay : Real;
    const by : Real;
    const N:Integer):Real;
var
    I   :   Integer;
    x   :   Real;
    Y   :   Real;
begin
    i:=1;
    Result:=0;
    repeat
        x:=ax+(bx-ax)*RandomReal();
        y:=ay+(by-ay)*RandomReal();
        if D(x,y) then
        begin
            Result:=Result+F(x,y);
        end;
        i:=i+1;
    until not(i<=n);
    Result:=(bx-ax)*(by-ay)*Result/N;
end;


end.

Бочканов Сергей, Быстрицкий Владимир
Copyright © 1999-2004
При поддержке проекта MANUAL.RU