Обчислення інтегралів методом Монте-Карло

Саратовського державного соціально-ЕКОНОМІЧНИЙ УНІВЕРСИТЕТ

ФАКУЛЬТЕТ ІНФОРМАТИКИ ТА ІНФОРМАЦІЙНИХ ТЕХНОЛОГІЙ

Курсова робота

Обчислення інтегралів МЕТОДОМ МОНТЕ - КАРЛО

Виконав:

Керівник:

Саратов, 2009

ЗМІСТ

ВСТУП

1. МАТЕМАТИЧНЕ ОБГРУНТУВАННЯ АЛГОРИТМУ ОБЧИСЛЕННЯ ІНТЕГРАЛА

1.1 Принцип роботи методу Монте - Карло

1.2 Застосування методу Монте - Карло для обчислення n - мірного інтеграла.

1.3 Сплайн - інтерполяція

1.4 Алгоритм розрахунку інтеграла

2. Генератор псевдовипадкових чисел

2.1 Генератор псевдовипадкових чисел стосовно методу Монте - Карло.

2.2 Алгоритм генератора псевдовипадкових чисел

2.3 Перевірка рівномірності розподілу генератора псевдовипадкових чисел.

ВИСНОВОК

СПИСОК ВИКОРИСТАНИХ ДЖЕРЕЛ

ВСТУП

Метою даної роботи є створення програмного продукту для участі в конкурсі, проведеному групою компаній «Траст» зі створення програмних розробок. Для реалізації було вибрано наступне технічної завдання:

Завдання 12 Обчислення інтегралів методом Монте - Карло.

Мета:

1. Реалізація генератора випадкових чисел для методу Монте - Карло.

2. Порівняння рівномірного розподілу і спеціально розробленого.

3. Обчислення тестового багатовимірного інтеграла в складній області.

Продукт:

1. Програмний код у вигляді функції на мові С + + або Fortran.

2. Тестові приклади у вигляді програми, що викликають реалізовані функції.

3. Огляд використаної літератури.

Для реалізації даного технічного завдання була обрана мова C + +. Код реалізований в інтегрованому середовищі розробки додатків Borland C + + Builder Enterprises і математично обгрунтований відповідний спосіб обчислення інтеграла.

1. МАТЕМАТИЧНЕ ОБГРУНТУВАННЯ АЛГОРИТМУ ОБЧИСЛЕННЯ ІНТЕГРАЛА

1.1 Принцип роботи методу Монте - Карло

Датою народження методу Монте - Карло визнано вважати 1949 рік, коли американські учені Н. Метрополіс і С. Услам опублікували статтю під назвою «Метод Монте - Карло», в якій були викладені принципи цього методу. Назва методу походить від назви міста Монте - Карло, який славився своїми гральними закладами, неодмінним атрибутом яких була рулетка - одне з найпростіших засобів отримання випадкових чисел з хорошим рівномірним розподілом, на використанні яких заснований цей метод.

Метод Монте - Карло це статистичний метод. Його використовують при обчисленні складних інтегралів, рішення систем алгебраїчних рівнянь високого порядку, моделюванні поведінки елементарних частинок, в теоріях передачі інформації, при дослідженні складних економічних систем.

Суть методу полягає в тому, що в завдання вводять випадкову величину , Що змінюється по якому те правилом . Випадкову величину вибирають таким чином, щоб бажана в задачі величина стала математичним очікування від , Тобто .

Таким чином, шукана величина визначається лише теоретично. Щоб знайти її чисельно необхідно скористатися статистичними методами. Тобто необхідно взяти вибірку випадкових чисел обсягом . Потім необхідно обчислити вибіркове середнє варіанту випадкової величини за формулою:

. (1)

Обчислене вибіркове середнє приймають за наближене значення .

Для отримання результату прийнятної точності необхідна велика кількість статистичних випробувань.

Теорія методу Монте - Карло вивчає способи вибору випадкових величин для вирішення різних завдань, а також способи зменшення дисперсії випадкових величин.

1.2 Застосування методу Монте - Карло для обчислення n - мірного інтеграла.

Розглянемо n - мірний інтеграл

для . (2)

Будемо вважати, що область інтегрування , І що обмежене безліч в . Отже, кожна точка х безлічі має n координат: .

Функцію візьмемо таку, що вона обмежена зверху і знизу на безлічі : .

Скористаємося обмеженістю безлічі і впишемо його в деякий n - мірний паралелепіпед , Наступним чином:

,

де - Мінімуми і максимуми, відповідно, - Ой координати всіх точок безлічі : .

Доопределяется подинтегральную функцію таким чином, щоб вона зверталася в нуль в точках паралелепіпеда , Які не належать :

(3)

Таким чином, рівняння (2) можна записати у вигляді

. (4)

Область інтегрування являє собою n - мірний паралелепіпед зі сторонами паралельними осям координат. Даний паралелепіпед можна однозначно задати двома вершинами , Які мають наймолодші та найстарші координати всіх точок паралелепіпеда.

Позначимо через n-мірний вектор, що має рівномірний розподіл у паралелепіпеді : , Де .

Тоді її щільність ймовірностей буде визначена таким чином

(5)

Значення подинтегральной функції від випадкового вектора буде випадковою величиною , Математичне очікування якої є середнім значенням функції на безлічі :

. (6)

Середнє значення функції на безлічі дорівнює відношенню значення шуканого інтеграла до обсягу паралелепіпеда :

(7)

Позначимо обсяг паралелепіпеда .

Таким чином, значення шуканого інтеграла можна виразити як добуток математичного сподівання функції і обсягу n - мірного паралелепіпеда :

(8)

Отже, необхідно знайти значення математичного очікування . Його наближене значення можна знайти провівши n випробувань, отримавши, таким чином, вибірку випадкових векторів, які мають рівномірний розподіл на . Позначимо і . Для оцінки математичного сподівання скористаємося результатом

, (9)

де ,

,

- Квантиль нормального розподілу, що відповідає довірчій ймовірності .

Помноживши подвійне нерівність з (9) на одержимо інтервал для I:

. (10)

Позначимо точкову оцінку . Отримуємо оцінку (з надійністю ):

. (11)

Аналогічно можна знайти вираз для відносної похибки :

. (12)

Якщо задана цільова абсолютна похибка , З (11) можна визначити обсяг вибірки, що забезпечує задану точність і надійність:

. (13)

Якщо задана цільова відносна похибка, з (12) отримуємо аналогічне вираз для об'єму вибірки:

. (14)

1.3 Сплайн - інтерполяція.

У даному програмному продукті реалізована можливість задавати додаткові обмеження області інтегрування двома двовимірними сплайн - поверхнями (для подинтегральной функції розмірності 3). Для завдання цих поверхонь використовуються двовимірні сплайни типу гнучкої пластинки \ 4 \.

Під сплайном (від англ. Spline - планка, рейка) зазвичай розуміють агрегатну функцію, збігається з функціями більш простої природи на кожному елементі розбиття своєї області визначення. Сплайн - функція має наступний вигляд:

. (15)

Вихідні дані є трійок точок .

Коефіцієнти і визначаються із системи:

, (16)

де ,

.

1.4 Алгоритм розрахунку інтеграла

Реалізований алгоритм включає наступні кроки:

1. вибирається початкове значення , Розігруються випадкові вектори з і визначаються і ;

2. в залежності від виду похибки (абсолютна, відносна) визначається досягнута похибка; якщо вона менше цільової, обчислення переривається;

3. за формулами (13) або (14) обчислюється новий обсяг вибірки;

4. обсяг вибірки збільшується на 20%

5. перехід до кроку 1;

6. кінець.

2. Генератор псевдовипадкових чисел

2.1 Генератор псевдовипадкових чисел стосовно методу Монте - Карло.

У будь-якому алгоритмі використовує метод Монте - Карло генератор псевдовипадкових чисел відіграє дуже важливу роль. Ступінь відповідності псевдовипадкових чисел заданому розподілу є важливим фактором проведення якісних статистичних випробувань.

2.2 Алгоритм генератора псевдовипадкових чисел

У програмі реалізований конгруентний метод генерації псевдовипадкових чисел \ 3 \:

, (17)

де = 8192,

= 67101323.

Авторський код, який реалізує захист від переповнення був, реалізований на С + +. Перед використання перші три числа послідовності видаляються. Для отримання чисел з інтервалу (0,1) усі числа діляться на .

2.3 Перевірка рівномірності розподілу генератора псевдовипадкових чисел.

Перевірка рівномірності розподілу псевдовипадкових чисел проводилася за допомогою стандартного критерію χ ² \ ² \.

Були використані 3 послідовності псевдовипадкових чисел, що визначаються стартовими значеннями 1, 1001, 1000000 довжиною 300000.

Інтервал (0,1) поділяються на 50 рівних інтервалів і програмно підраховувалися абсолютні частоти (рис. 1).

Рис. 1

Результати перевірки наведено в Таблиці 1.

Таблиця 1

Отже, рівномірність розподілу не відкидається на рівні 5%.

ВИСНОВОК

На закінчення можна сказати, що поставлене завдання було повністю виконано. Тобто на мові С + + було розроблено генератор псевдовипадкових чисел, функція розраховує інтеграл методом Монте - Карло (Додаток 1); був проведений розрахунок тестових багатовимірних інтегралів (Додаток 2); в інтегрованому середовищі розробки додатків Borland C + + Builder Enterprises 7.0 був створений програмний продукт «CarloS», який реалізує описані вище алгоритми (Додаток 3).

СПИСОК ВИКОРИСТАНИХ ДЖЕРЕЛ

1. Бережна Є. В., Бережний В. І. Математичні методи моделювання економічних систем. - М.: Фінанси і статистика, 2001. - 368 с.

2. Мюллер П., Нойман П., Шторм Р. Таблиці з математичної статистики. - М.: Фінанси і статистика, 1982. - 278 с.

3. Теннант-Сміт Дж. Бейсік для статистиків. - М.: Мир, 1988. - 208 с.

4. Baranger J. Analyse numérique. Hermann, 1991.

5. Маделунга Е. Математичний апарат фізики. Довідкове керівництво. М.: Наука, 1968., С.287.

6. В.Є. Гмурман Теорія ймовірностей і математична статистика - М.: Вища школа, 2003

ДОДАТОК 1

Лістинг ОСНОВНИХ ФУНКЦІЙ

Лістинг 1 Функція розрахунку інтеграла

void integral ()

{

/ / Обчислення інтеграла методом Монте - Карло

/ / Розмірність області інтегрування

unsigned d_int = fun_dim;

//----- 3 d графік ---------------------------------------- ----------------

/ / Максимальне число трійок

unsigned plot_dim_max = 10000;

/ / Матриця трійок

pmatd xyz, xyz_tmp;

if (d_int == 3) xyz = new matd (plot_dim_max, 3);

//------------------------------------------------ -------------------------

/ / Індикатор відносної похибки

mcres.relok = Read1double ("error_type.txt");

/ / Цільова похибка

mcres.dlt_int = Read1double ("error_value.txt");

/ / Номер стандартного значення довірчої ймовірності (починаючи з 0)

int nome_int = Read1double ("error_omega.txt");

/ / ГВЧ

unsigned long b = m_rng * m_rng-d_rng, c, r, i, PSChunk;

/ / "Паросток" ГВЧ

mcres.rng_seed = Read1double ("rng_seed.txt");

pmatd fun_b, fun_A, con_b, con_A, con_U, con_v, \

a_int, b_int, ba_int, x_int, xyz_top, xyz_bottom;

unsigned j, ii, jj, con_ok;

struct date dat;

struct time tim;

pspl2d sp_top, sp_bottom;

/ / Квантилі нормального розподілу

double omegas _ int [6] = {0.9,0.95,0.99,0.999,0.9999,0.99999};

double zs_int [6] = {1.64485362695147,1.95996398454005,2.5758293035489, \

3.29052673149191, 3.89059188641317, 4.4171734134667};

mcres.omega_int = omegas_int [nome_int];

mcres.z_int = zs_int [nome_int];

double fun_cd, con_wd, fu_int, con_sum, sum1_int, sum2_int;

/ / Вид інтегрованої функції

/ / 0 - постійна

/ / 1 - лінійна

/ / 2 - квадратична

mcres.fun_type = Read1double ("fun_kind.txt");

/ / Вид системи обмежень

/ / 0 - відсутні (весь паралелепіпед)

/ / 1 - лінійні

/ / 2 - квадратичне

/ / 3 - сплайн - поверхні

mcres.con_type = Read1double ("con_type.txt");

/ / Завантаження параметрів інтегрованої функції

switch (mcres.fun_type)

{

case 2: fun_A = new matd ("fun_A.txt");

case 1: fun_b = new matd ("fun_b.txt");

case 0: fun_cd = Read1double ("fun_c.txt");

}

/ / Завантаження параметрів обмежень

switch (mcres.con_type)

{

case 3: / / сплайн - поверхні

/ / Верхня

xyz_top = new matd ("xyz_top.txt");

/ / Нижня

xyz_bottom = new matd ("xyz_bottom.txt");

/ / Двовимірна інтерполяція

sp_top = new spl2d (xyz_top);

sp_bottom = new spl2d (xyz_bottom);

break;

case 2: / / квадратична функція обмежень

con_U = new matd ("con_U.txt");

con_v = new matd ("con_v.txt");

con_wd = Read1double ("con_w.txt");

break;

case 1: / / лінійні обмеження

con_b = new matd ("con_b.txt"); con_A = new matd ("con_A.txt");

}

/ / Осяжний паралелепіпед

a_int = new matd ("con_xmin.txt");

b_int = new matd ("con_xmax.txt");

/ / Різниця меж паралелепіпеда

ba_int = new matd;

ba_int = & (* b_int - (* a_int));

/ / Аргумент інтегрованої функції

x_int = new matd (d_int, 1);

/ / Обсяг осяжний паралелепіпеда

mcres.V0_int = 1;

for (j = 1; j <= d_int; j + +)

{

if (_p (ba_int, j, 1) <= 0)

{

DbBox ("Нижня межа осяжний паралелепіпеда вище верхньої для \

координати ", j);

goto clean_exit;

}

mcres.V0_int = mcres.V0_int * _p (ba_int, j, 1);

}

/ / Початковий обсяг вибірки

mcres.n1_int = 10000;

/ / Основний цикл для досягнення заданої точності

/ / Число ітерацій, що необхідний для досягнення заданої точності

mcres.n_ite = 0;

getdate (& dat); gettime (& tim); mcres.t_start = dostounix (& dat, & tim);

WaitForm-> Show ();

while (1)

{

mcres.n_ite + +;

WaitForm-> Edit1-> Text = mcres.n_ite;

WaitForm-> Edit2-> Text = mcres.n1_int;

WaitForm-> ProgressBar1-> Position = 0;

WaitForm-> Refresh ();

/ / Генерація випадкових точок і накопичення суми

sum1_int = 0; sum2_int = 0;

mcres.in_G_int = 0;

PSChunk = long (mcres.n1_int/50.0);

/ / Запуск ГВЧ

r = mcres.rng_seed;

for (i = 1; i <3; i + +)

{

c = int (r / m_rng);

r = b * c + m_rng * (r-m_rng * c);

if (r> d_rng) r = r-d_rng;

}

for (i = 1; i <= mcres.n1_int; i + +)

{

/ / Випадковий вектор

for (j = 1; j <= d_int; j + +)

{

/ / Випадкове число

c = int (r / m_rng);

r = b * c + m_rng * (r-m_rng * c);

if (r> d_rng) r = r-d_rng;

_p (x_int, j, 1) = _p (a_int, j, 1) + _p (ba_int, j, 1) * double (r) / d_rng;

}

/ / Прогрес

if (! (i% PSChunk))

{

WaitForm-> ProgressBar1-> Position = 100.0 * (i-1) / (mcres.n1_int-1);

WaitForm-> Refresh ();

}

/ / Перевірка обмеження

con_ok = 1;

switch (mcres.con_type)

{

case 3: / / сплайн - поверхні

if ((_p (x_int, 3,1) <sp_bottom-> f (_p (x_int, 1,1), \

_p (x_int, 2,1 )))||(_ p (x_int, 3,1)> sp_top-> f (_p (x_int, 1,1), _p (x_int, 2,1)))) con_ok = 0 ;

break;

case 2: / / квадратична функція обмежень

con_sum = 0;

for (ii = 1; ii <= d_int; ii + +)

for (jj = 1; jj <= d_int; jj + +)

if (_p (con_U, ii, jj)! = 0)

con_sum + = _p (x_int, ii, 1) * _p (con_U, ii, jj) * _p (x_int, jj, 1);

for (ii = 1; ii <= d_int; ii + +)

if (_p (con_v, ii, 1)! = 0)

con_sum + = _p (con_v, ii, 1) * _p (x_int, ii, 1);

if (con_sum> con_wd) con_ok = 0;

break;

case 1: / / лінійна функція обмежень

for (ii = 1; ii <= con_A-> nl; ii + +)

{

con_sum = 0;

for (jj = 1; jj <= d_int; jj + +)

con_sum + = _p (con_A, ii, jj) * _p (x_int, jj, 1);

if (con_sum> _p (con_b, ii, 1)) {con_ok = 0; break;}

}

}

fu_int = 0;

if (con_ok! = 0)

{

mcres.in_G_int + +;

/ / Точки 3d графіка

if (d_int == 3)

if (mcres.in_G_int <= plot_dim_max)

{

_p (xyz, mcres.in_G_int, 1) = _p (x_int, 1,1);

_p (xyz, mcres.in_G_int, 2) = _p (x_int, 2,1);

_p (xyz, mcres.in_G_int, 3) = _p (x_int, 3,1);

}

/ / Значення інтегрованої функції

switch (mcres.fun_type)

{

case 2: / / квадратичний член

for (ii = 1; ii <= d_int; ii + +)

for (jj = 1; jj <= d_int; jj + +)

if (_p (fun_A, ii, jj)! = 0)

fu_int + = _p (x_int, ii, 1) * _p (fun_A, ii, jj) * _p (x_int, jj, 1);

case 1: / / лінійний член

for (ii = 1; ii <= d_int; ii + +)

if (_p (fun_b, ii, 1)! = 0)

fu_int + = _p (fun_b, ii, 1) * _p (x_int, ii, 1);

case 0: / / постійна

fu_int + = fun_cd;

}

}

sum1_int + = fu_int; sum2_int + = fu_int * fu_int;

}

/ / Оцінка мат. очікування і дисперсії

mcres.f1_int = sum1_int/mcres.n1_int;

mcres.vari_int = (sum2_int-sum1_int * sum1_int/mcres.n1_int) / (mcres.n1_int-1);

/ / Розрахунок похибки

if (mcres.relok == 0)

{

/ / Абсолютна похибка

mcres.deltar = mcres.V0_int * mcres.z_int * sqrt (mcres.vari_int/mcres.n1_int);

}

else

{

/ / Відносна похибка

if (mcres.f1_int! = 0)

{

mcres.deltar = mcres.z_int / fabs (mcres.f1_int) * sqrt (mcres.vari_int/mcres.n1_int);

}

else

{

/ / Форма результатів

mcres. inte _ int = 0;

mcres.deltar = 0;

getdate (& dat); gettime (& tim); mcres.t_end = dostounix (& dat, & tim);

mcres.t_calc = mcres.t_end-mcres.t_start;

InfoBox ("Оцінка інтеграла = 0 (обрана относ. Похибка), обчислення \

перервано. ");

ResultForm-> Show ();

WaitForm-> Close ();

goto clean_exit;

}

}

WaitForm-> Edit3-> Text = mcres.deltar;

WaitForm-> Refresh ();

if (mcres.deltar <mcres.dlt_int)

{

/ / Точність достатня

mcres.inte_int = mcres.V0_int * mcres.f1_int;

getdate (& dat); gettime (& tim); mcres.t_end = dostounix (& dat, & tim);

mcres.t_calc = mcres.t_end-mcres.t_start;

ResultForm-> Show ();

break;

}

/ / Обчислення нового обсягу вибірки

if (mcres.relok == 0)

{

/ / Абс. Похибка

mcres.n1_int = ceil (mcres.vari_int * pow (mcres.V0_int * mcres.z_int / mcres.dlt_int, 2));

}

else

{

/ / Отн. Похибка

mcres.n1_int = ceil (mcres.vari_int * pow (mcres.z_int/mcres.dlt_int/mcres.f1_int, 2));

}

/ / Коректування обсягу вибірки в більшу сторону

/ / Для скорочення числа ітерацій

mcres.n1_int = 1.2 * mcres.n1_int;

/ / Мінімальний обсяг вибірки

if (mcres.n1_int <1000) mcres.n1_int = 1000;

} / / Кінець основного циклу

WaitForm-> Close ();

/ / 3d графік

if (d_int == 3)

{

if (mcres.in_G_int == 0)

{

/ / Безліч точок порожньо

Zero_File ("xyz.txt");

}

else

if (mcres.in_G_int <xyz-> nl)

{

/ / Точок не набралося, щоб заповнити матрицю

xyz_tmp = new matd (mcres.in_G_int, 3);

for (i = 1; i <= mcres.in_G_int; i + +)

{

_p (xyz_tmp, i, 1) = _p (xyz, i, 1);

_p (xyz_tmp, i, 2) = _p (xyz, i, 2);

_p (xyz_tmp, i, 3) = _p (xyz, i, 3);

}

xyz_tmp-> txprint ("xyz.txt");

delete xyz_tmp;

}

else

{

/ / Вся матриця заповнена

xyz-> txprint ("xyz.txt");

}

} / / Кінець d_int == 3

clean_exit:

/ / Очищення пам'яті

if (d_int == 3) delete xyz;

switch (mcres.fun_type)

{

case 2: delete fun_A;

case 1: delete fun_b;

}

switch (mcres.con_type)

{

case 3: delete xyz_top, xyz_bottom, sp_top, sp_bottom; break;

case 2: delete con_U, con_v; break;

case 1: delete con_b, con_A;

}

delete a_int, b_int, ba_int, x_int;

} / / Integral

Листинг 2 структура для зберігання результатів розрахунку інтеграла

struct mcres_struct

{

/ / Індикатор відносної похибки

int relok;

/ / Цільова похибка

double dlt_int;

/ / Досягнута похибка

double deltar;

/ / Довірча ймовірність

double omega_int;

/ / Квантиль норм. розподілу

double z_int;

/ / "Паросток" ГВЧ

unsigned long rng_seed;

/ / ÷ число ітерацій, що необхідний для досягнення заданої точності

unsigned n_ite;

/ / Обсяг вибірки на останній ітерації

unsigned long n1_int;

/ / Число точок потрапили в область інтегрування

unsigned in_G_int;

/ / Інтеграл

double inte_int;

/ / Обсяг осяжний паралелепіпеда

double V0_int;

/ / Вибіркове середнє

double f1_int;

/ / Вибіркова дисперсія

double vari_int;

/ / Час початку рахунки

time_t t_start;

/ / Час закінчення рахунку

time_t t_end;

/ / Тривалість обчислення інтеграла

time_t t_calc;

/ / Вид інтегрованої функції

int fun_type;

/ / Вид системи огрніченій

int con_type;

}; / / Mcres_struct

ДОДАТОК 2

ТЕСТОВІ ПРИКЛАДИ

Приклад 1 Інтеграл від квадратичної функції по 3-мірному симплекс.

Точне значення інтеграла:

Наближене значення знайдено для цільової абсолютної похибки 0.00001.

Похибка: 0.000034416630896 або 0.014749984670%.

Приклади 2-10 Обсяги багатовимірних куль

Точні і наближені обсяги багатовимірних куль наведено в таблиці нижче.