Студопедия

КАТЕГОРИИ:

АстрономияБиологияГеографияДругие языкиДругоеИнформатикаИсторияКультураЛитератураЛогикаМатематикаМедицинаМеханикаОбразованиеОхрана трудаПедагогикаПолитикаПравоПсихологияРиторикаСоциологияСпортСтроительствоТехнологияФизикаФилософияФинансыХимияЧерчениеЭкологияЭкономикаЭлектроника


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




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

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

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

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

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

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

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

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

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

. (9.9)

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

 

(9.10)

 

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

 

, (9.11)

где xi равномерно распределены на отрезке [a, b], yi равномерно распределены на отрезке [c, d], а zi равномерно распределены на отрезке [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, а в столбце BN = 1000000, поэтому во втором столбце получены более точные значения.

 

Таблица 9.3.

  A B
4,513159 4,998266
4,867957 5,002396
5,523045 4,998898
4,575598 5,001704
4,303717 5,003065

 

Составим программу на 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 испытаний для интеграла

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

 


Поделиться:

Дата добавления: 2015-09-13; просмотров: 159; Мы поможем в написании вашей работы!; Нарушение авторских прав





lektsii.com - Лекции.Ком - 2014-2024 год. (0.007 сек.) Все материалы представленные на сайте исключительно с целью ознакомления читателями и не преследуют коммерческих целей или нарушение авторских прав
Главная страница Случайная страница Контакты