САРАТОВСЬКИЙ ДЕРЖАВНИЙ СОЦІАЛЬНО-ЕКОНОМІЧНИЙ УНІВЕРСИТЕТ
ФАКУЛЬТЕТ ІНФОРМАТИКИ ТА ІНФОРМАЦІЙНИХ ТЕХНОЛОГІЙ
Курсова робота
Обчислення інтегралів МЕТОДОМ МОНТЕ - КАРЛО
Виконав:
Керівник:
Саратов, 2009
ЗМІСТ
ВСТУП
1. МАТЕМАТИЧНЕ ОБГРУНТУВАННЯ АЛГОРИТМУ ОБЧИСЛЕННЯ ІНТЕГРАЛА
1.1 Принцип роботи методу Монте - Карло
1.2 Застосування методу Монте - Карло для обчислення n - мірного інтеграла.
1.3 Сплайн - інтерполяція 8
1.4 Алгоритм розрахунку інтеграла
2. Генератор псевдовипадкових чисел
2.1 Генератор псевдовипадкових чисел стосовно до методу Монте - Карло.
2.2 Алгоритм генератора псевдовипадкових чисел
2.3 Перевірка рівномірності розподілу генератора псевдовипадкових чисел.
ВИСНОВОК
СПИСОК ВИКОРИСТАНИХ ДЖЕРЕЛ
Метою даної роботи є створення програмного продукту для участі в конкурсі, проведеному групою компаній «Траст» зі створення програмних розробок. Для реалізації було вибрано наступне технічної завдання:
Завдання 12 Обчислення інтегралів методом Монте - Карло.
Мета:
1) Реалізація генератора випадкових чисел для методу Монте - Карло.
2) Порівняння рівномірного розподілу і спеціально розробленого.
3) Обчислення тестового багатовимірного інтеграла в складній області.
Продукт:
1) Програмний код у вигляді функції на мові С + + або Fortran.
2) Тестові приклади у вигляді програми, що викликають реалізовані функції.
3) Огляд використаної літератури.
Для реалізації даного технічного завдання була обрана мова C + +. Код реалізований в інтегрованому середовищі розробки додатків Borland C + + Builder Enterprises і математично обгрунтовано відповідний спосіб обчислення інтеграла.
Метод Монте - Карло це статистичний метод. Його використовують при обчисленні складних інтегралів, рішення систем алгебраїчних рівнянь високого порядку, моделюванні поведінки елементарних частинок, в теоріях передачі інформації, при дослідженні складних економічних систем.
Суть методу полягає в тому, що в завдання вводять випадкову величину , Що змінюється по якому те правилом . Випадкову величину вибирають таким чином, щоб бажана в задачі величина стала математичним очікування від , Тобто .
Таким чином, шукана величина визначається лише теоретично. Щоб знайти її чисельно необхідно скористатися статистичними методами. Тобто необхідно взяти вибірку випадкових чисел об'ємом . Потім необхідно обчислити вибіркове середнє варіанту випадкової величини за формулою:
. (1)
Обчислення вибіркове середнє приймають за наближене значення .
Для отримання результату прийнятної точності необхідна велика кількість статистичних випробувань.
Теорія методу Монте - Карло вивчає способи вибору випадкових величин для вирішення різних завдань, а також способи зменшення дисперсії випадкових величин.
для . (2)
Будемо вважати, що область інтегрування , І що обмежене безліч в . Отже, кожна точка х безлічі має n координат: .
Функцію візьмемо таку, що вона обмежена зверху і знизу на безлічі : .
Скористаємося обмеженістю безлічі і впишемо його в деякий n - мірний паралелепіпед , Наступним чином:
,
де - Мінімуми і максимуми, відповідно, - Ой координати всіх точок множини : .
Доопределяют підінтегральна функція таким чином, щоб вона зверталася в нуль в точках паралелепіпеда , Які не належать :
(3)
Таким чином, рівняння (2) можна записати у вигляді
. (4)
Область інтегрування являє собою n - мірний паралелепіпед зі сторонами паралельними осям координат. Даний паралелепіпед можна однозначно задати двома вершинами , Які мають наймолодші та найстарші координати всіх точок паралелепіпеда.
Позначимо через n-мірний вектор, що має рівномірний розподіл у паралелепіпеді : , Де .
Тоді її щільність ймовірностей буде визначена наступним чином
(5)
Значення підінтегральної функції від випадкового вектора буде випадковою величиною , Математичне сподівання якої є середнім значенням функції на множині :
. (6)
Середнє значення функції на множині дорівнює відношенню значення шуканого інтеграла до обсягу паралелепіпеда :
(7)
Позначимо обсяг паралелепіпеда .
Таким чином, значення шуканого інтеграла можна виразити як добуток математичного сподівання функції та обсягу n - мірного паралелепіпеда :
(8)
Отже, необхідно знайти значення математичного очікування . Його наближене значення можна знайти провівши n випробувань, отримавши, таким чином, вибірку випадкових векторів, що мають рівномірний розподіл на . Позначимо і . Для оцінки математичного сподівання скористаємося результатом
, (9)
де ,
,
- Квантиль нормального розподілу, що відповідає довірчій ймовірності .
Помноживши подвійне нерівність з (9) на матимемо інтервал для I:
. (10)
Позначимо точкову оцінку . Отримуємо оцінку (з надійністю ):
. (11)
Аналогічно можна знайти вираз для відносної похибки :
. (12)
Якщо задана цільова абсолютна похибка , З (11) можна визначити обсяг вибірки, що забезпечує задану точність і надійність:
. (13)
Якщо задана цільова відносна похибка, з (12) отримуємо аналогічне вираз для обсягу вибірки:
. (14)
Під сплайном (від англ. Spline - планка, рейка) зазвичай розуміють агрегатну функцію, збігається з функціями більш простої природи на кожному елементі розбиття своїй області визначення. Сплайн - функція має наступний вигляд:
. (15)
Вихідні дані представляють собою трійок точок .
Коефіцієнти і визначаються з системи:
, (16)
де ,
.
1) вибирається початкове значення , Розігруються випадкові вектори з і визначаються і ;
2) залежно від виду похибки (абсолютна, відносна) визначається досягнута похибка; якщо вона менше цільової, обчислення переривається;
3) за формулами (13) або (14) обчислюється новий обсяг вибірки;
4) обсяг вибірки збільшується на 20%
5) перехід до кроку 1;
6) кінець.
, (17)
де = 8192,
= 67101323.
Авторський код, що реалізовує захист від переповнення був, реалізований на С + +. Перед використання перші три числа послідовності видаляються. Для отримання чисел з інтервалу (0,1) всі числа діляться на .
Були використані 3 послідовності псевдовипадкових чисел, що визначаються стартовими значеннями 1, 1001, 1000000 довжиною 300000.
Інтервал (0,1) подразделялся на 50 рівних інтервалів і програмно підраховувалися абсолютні частоти (рис. 1).
Рис. 1
Результати перевірки наведено в Таблиці 1.
Таблиця 1
Отже, рівномірність розподілу не відкидається на рівні 5%.
2. Мюллер П., Нойман П., Шторм Р. Таблиці з математичної статистики. - М.: Фінанси і статистика, 1982. - 278 с.
3. Теннант-Сміт Дж. Бейсік для статистиків. - М.: Мир, 1988. - 208 с.
4. Baranger J. Analyse numérique. Hermann, 1991.
5. Маделунга Е. Математичний апарат фізики. Довідкове керівництво. М.: Наука, 1968., С.287.
6. В.Є. Гмурман Теорія ймовірностей і математична статистика - М.: Вища школа, 2003
ФАКУЛЬТЕТ ІНФОРМАТИКИ ТА ІНФОРМАЦІЙНИХ ТЕХНОЛОГІЙ
Курсова робота
Обчислення інтегралів МЕТОДОМ МОНТЕ - КАРЛО
Виконав:
Керівник:
Саратов, 2009
ЗМІСТ
ВСТУП
1. МАТЕМАТИЧНЕ ОБГРУНТУВАННЯ АЛГОРИТМУ ОБЧИСЛЕННЯ ІНТЕГРАЛА
1.1 Принцип роботи методу Монте - Карло
1.2 Застосування методу Монте - Карло для обчислення n - мірного інтеграла.
1.3 Сплайн - інтерполяція 8
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.2 Застосування методу Монте - Карло для обчислення n - мірного інтеграла.
Розглянемо n - мірний інтегралБудемо вважати, що область інтегрування
Функцію
Скористаємося обмеженістю безлічі
де
Доопределяют підінтегральна функція
Таким чином, рівняння (2) можна записати у вигляді
Область інтегрування являє собою n - мірний паралелепіпед
Позначимо через
Тоді її щільність ймовірностей
Значення підінтегральної функції
Середнє значення функції на множині
Позначимо
Таким чином, значення шуканого інтеграла можна виразити як добуток математичного сподівання функції та обсягу n - мірного паралелепіпеда
Отже, необхідно знайти значення математичного очікування
де
Помноживши подвійне нерівність з (9) на
Позначимо
Аналогічно можна знайти вираз для відносної похибки
Якщо задана цільова абсолютна похибка
Якщо задана цільова відносна похибка, з (12) отримуємо аналогічне вираз для обсягу вибірки:
1.3 Сплайн - інтерполяція.
У даному програмному продукті реалізована можливість задавати додаткові обмеження області інтегрування двома двовимірними сплайн - поверхнями (для підінтегральної функції розмірності 3). Для завдання цих поверхонь використовуються двовимірні сплайни типу гнучкої пластинки \ 4 \.Під сплайном (від англ. Spline - планка, рейка) зазвичай розуміють агрегатну функцію, збігається з функціями більш простої природи на кожному елементі розбиття своїй області визначення. Сплайн - функція має наступний вигляд:
Вихідні дані представляють собою
Коефіцієнти
де
1.4 Алгоритм розрахунку інтеграла
Реалізований алгоритм включає наступні кроки:1) вибирається початкове значення
2) залежно від виду похибки (абсолютна, відносна) визначається досягнута похибка; якщо вона менше цільової, обчислення переривається;
3) за формулами (13) або (14) обчислюється новий обсяг вибірки;
4) обсяг вибірки збільшується на 20%
5) перехід до кроку 1;
6) кінець.
2. Генератор псевдовипадкових чисел
2.1 Генератор псевдовипадкових чисел стосовно до методу Монте - Карло.
У будь-якому алгоритмі використовує метод Монте - Карло генератор псевдовипадкових чисел грає дуже важливу роль. Ступінь відповідності псевдовипадкових чисел заданого розподілу є важливим чинником проведення якісних статистичних випробувань.2.2 Алгоритм генератора псевдовипадкових чисел
У програмі реалізований конгруентний метод генерації псевдовипадкових чисел \ 3 \:де
Авторський код, що реалізовує захист від переповнення був, реалізований на С + +. Перед використання перші три числа послідовності видаляються. Для отримання чисел з інтервалу (0,1) всі числа діляться на
2.3 Перевірка рівномірності розподілу генератора псевдовипадкових чисел.
Перевірка рівномірності розподілу псевдовипадкових чисел проводилася за допомогою стандартного критерію χ 2 \ 2 \.Були використані 3 послідовності псевдовипадкових чисел, що визначаються стартовими значеннями 1, 1001, 1000000 довжиною 300000.
Інтервал (0,1) подразделялся на 50 рівних інтервалів і програмно підраховувалися абсолютні частоти (рис. 1).
Рис. 1
Результати перевірки наведено в Таблиці 1.
Таблиця 1
стартове значення ГВЧ | |||
1 | 1001 | 1000000 | |
хі-квадрат | 44.0533333333333 | 45.007 | 48.618 |
df | 50 | 50 | 50 |
p-значення | 0.709735881642893 | 0.673522612551685 | 0.528941919633451 |
ВИСНОВОК
На закінчення можна сказати, що поставлене завдання було повністю виконано. Тобто на мові С + + було розроблено генератор псевдовипадкових чисел, функція розраховує інтеграл методом Монте - Карло (Додаток 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;
Лістинг ОСНОВНИХ ФУНКЦІЙ
Лістинг 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 Обсяги багатовимірних куль
Точні і наближені обсяги багатовимірних куль приведені в наступній таблиці.
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.000034416630896 або 0.014749984670%.
Приклади 2-10 Обсяги багатовимірних куль
Точні і наближені обсяги багатовимірних куль
Обсяг точний [1] | Обсяг | Оцінка CarloS [3] | Відносна похибка,% | |
2 | 3.1415926535897932385 | 3.1504 | 0.280346543342 | |
3 | 4.1887902047863909846 | 4.2032 | 0.344008520578 | |
4 | 4.9348022005446793096 | 4.98099547511312 | .936071451118 | |
5 | 5.2637890139143245968 | 5.18913116403891 | -1.4183290720439 | |
6 | 5.1677127800499700296 | 5.16153372226575 | -. 1195704569352 | |
7 | 4.7247659703314011698 | 4.70163814726423 | -. 4895019819476 | |
8 | 4.0587121264167682184 | 3.98117943332154 | -1.9102782035357 | |
9 | 3.2985089027387068695 | 3.30542485033746 | .209668908064 | |
10 | 2.5501640398773454440 | 2.55096385956571 | .31363460384 E-1 |
[1] Джерело [5], с. 287.
[2] Обчислено в Maple (20 значущих цифр).
[3] Розрахунок з цільовою відносною похибкою 2%