Вычисление кратных интегралов

Метод Монте-Карло для вычисления одномерных интегралов обычно не применяется, так как для получения высокой точности более удобны квадратурные формулы. Этот метод оказывается более эффективным при вычислении кратных интегралов, когда кубатурные формулы для получения достижения малой погрешности слишком громоздки и требуют большого объема вычислений.

При использовании квадратурных или кубатурных формул, число операций быстро возрастает с ростом размерности интеграла. Например, если для вычисления одномерного интеграла методом трапеций с заданной точностью необходимо вычислить сумму порядка N слагаемых, то для вычисления двойного интеграла тем же методом необходимо сложить порядка N² слагаемых, а для тройного интеграла число слагаемых составляет порядка N³.

Число испытаний N, требующихся для достижения заданной точности ε приближенного значения, в методе Монте-Карло есть величина порядка и не зависит от размерности интеграла.

Приведем критерий выбора между кубатурной формулой p-го порядка точности и методом Монте-Карло для вычисления с точностью ε кратного интеграла функции m переменных (обоснование содержится в [9]):

Если число измерений m < 2p, то выгодней использовать кубатурные или квадратурные формулы, если m > 2p — лучше метод Монте-Карло.

Например, если p = 1, тройные интегралы выгоднее вычислять методом Монте-Карло, а одномерные — квадратурными формулами.

Если p = 2, то выгоднее вычислять методом Монте-Карло пятимерные интегралы, а одномерные, двойные и тройные — квадратурными или кубатурными формулами.

Рассмотрим конкретные формулы метода Монте-Карло для вычисления кратных интегралов, получающиеся способом, который применялся для вывода формулы (9.7).

Пусть требуется вычислить двойной интеграл

. (9.9)

Проведем серию из N испытаний случайной точки , где x_i равномерно распределены на отрезке [a, b], а y_i равномерно распределены на отрезке [c, d]. Вычислим интеграл (9.9) по формуле

(9.10)

Для тройного интеграла аналогично получим формулу

, (9.11)

где x_i равномерно распределены на отрезке [a, b], y_i равномерно распределены на отрезке [c, d], а z_i равномерно распределены на отрезке [p, q]; N — число испытаний.

Для m-кратного интеграла формула метода Монте-Карло имеет вид

(9.12)

Пример 9.2. Вычислим двойной интеграл

.

Решение в программе Excel. Точное значение равно 5.

Построим макрос — функцию для вычисления интеграла по формуле (9.10). С помощью команды меню «Сервис — Макросы — Редактор Visual Basic» откроем окно редактора, выполним команду «Insert — Module» и введем тексты описаний функций:

Function f(ByVal x, ByVal y)

f = x * y

End Function

Function int_MK2(ByVal a, ByVal b, ByVal c, ByVal d, ByVal N)

Randomize

s = 0

For i = 1 To N

x = Rnd * (b - a) + a

y = Rnd * (d - c) + c

s = s + f(x, y)

Next i

int_MK2 = (b - a) * (d - c) * s / N

End Function

Введем в ячейку A1 формулу =int_MK2(0;2;2;3;100), а в ячейку B1 — =int_MK2(0;2;2;3;1000000).

Выделим диапазон A1:B1 и протянем маркер заполнения до строки 5. Получим таблицу 9.3.

В столбце A значения соответствуют сериям испытаний с объемом N = 100, а в столбце B — N = 1000000, поэтому во втором столбце получены более точные значения.

Таблица 9.3.

Составим программу на C++ для вычисления тройного интеграла по формуле (9.12).

#include <iostream.h>

#include <stdlib.h>

#include <math.h>

double f(double x, double y, double z);

typedef double (*PF)(double,double,double);

double MK3(PF f,double a1, double b1, double a2, double b2,

double a3, double b3, const int n);

int main(){

double a1, b1, a2, b2,a3, b3, y; PF pf; int n;

cout << " a1 = "; cin >> a1; cout << " b1 = "; cin >> b1;

cout << " a2 = "; cin >> a2; cout << " b2 = "; cin >> b2;

cout << " a3 = "; cin >> a3; cout << " b3 = "; cin >> b3;

cout << " n = "; cin >> n;

pf = f;

y = MK3(pf, a1, b1, a2, b2, a3, b3, n);

cout << " Integral = " << y;

cout << "\n Press any key & Enter "; cin >> n; // pause

return 0;

}

double f(double x, double y, double z){

double r;

r = x*y*z;

return r ;

}

double MK3(PF f,double a1, double b1, double a2, double b2,

double a3, double b3, const int n){

double s,x1,x2,x3; int i;

s = 0; randomize();

for(i=1; i<=n; i++){

x1 = a1+(b1-a1)*random(100)/100; x2 = a2+(b2-a3)*random(100)/100;

x3 = a3+(b3-a3)*random(100)/100; s = s+f(x1,x2,x3);

}

s = (b1-a1)*(b2-a2)*(b3-a3)*s/n;

return s;

}

Результат при n = 1000000 испытаний для интеграла

приведен ниже: