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

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

ВВЕДЕНИЕ

Целью данной работы является создание программного продукта для участия в конкурсе, проводимом группой компаний «Траст» по созданию программных разработок. Для реализации было выбрано следующее технической задание:Задание Вычисление интегралов методом Монте — Карло.
Цель
:
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) можно записать в виде .

7 стр., 3046 слов

Доклад: Применение экономико-математических методов в экономике

... века ознаменовались широким использованием математики в экономике. В XX в. математические методы моделирования используются столь широко, ... программирование, сетевые методы, метод Монте-Карло (метод статистических испытаний), методы теории надежности, случайные процессы, марковские ... экономики или поведение одной из составляющих в среде остальных. Основные объекты приложения моделирования в ...

(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 Проверка равномерности распределения генератора псевдосл
учайных чисел. Проверка равномерности распределения псевдослучайных чисел проводилась с помощью стандартного критерия ч2
\2\.Были использованы 3 последовательности псевдослучайных чисел, определяемых стартовыми значениями 1, 1001, 1000000 длиной 300000.Интервал (0,1) подразделялся на равных интервалов и программно подсчитывались абсолютные частоты (рис. 1).

Рис. 1

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

Таблица 1

стартовое значение ГСЧ

1

1001

1000000

хи-квадрат

44.0533333333333

45.007

48.618

p-значение

0.709735881642893

0.673522612551685

0.528941919633451

Следовательно, равномерность распределения не отвергается на уровне 5%.

ЗАКЛЮЧЕНИЕ


В заключение можно сказать, что поставленная задача была полностью выполнена. То есть на языке С++ были разработаны генератор псевдослучайных чисел, функция рассчитывающая интеграл методом Монте — Карло (Приложение 1); был проведен расчет тестовых многомерных интегралов (Приложение 2); в интегрированной среде разработки приложений Borland C++ Builder Enterprises 7.0 был создан программный продукт «CarloS», реализующий описанные выше алгоритмы (Приложение 3).СПИСОК ИСПОЛЬЗОВАННЫХ ИСТОЧНИКОВ

1. Бережная Е. В., Бережной В. И. Математические методы моделирования экономических систем. — М.: Финансы и статистика, 2001. — 368 с.

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

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

4. Baranger J. Analyse numerique. Hermann, 1991.

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

6. В.Е. Гмурман Теория вероятностей и математическая статистика — М.: Высшая школа, 2003

1. ПРИЛОЖЕНИЕ 1

ЛИСТИНГИ ОСНОВНЫХ ФУНКЦИЙ

Листинг 1 Функция расчета интеграла


void

integral ()


{


вычисление интеграла методом Монте

Карло


размерность области интегрирования


unsigned

d_int=fun_dim;


//—— 3 d

график


/i> максимальное число троек


unsigned

plot_dim_max=10000;


/i> матрица троек

pmatd xyz,xyz_tmp;

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


//————————————————————————-


/i> индикатор относительной погрешности

mcres.relok=Read1double(«error_type.txt»);


/i> целевая погрешность

mcres.dlt_int=Read1double(«error_value.txt»);


номер стандартного значения доверительной вероятности (нач

и
ная с 0)


int

nome_int=Read1double(«error_omega.txt»);


/i> ГСЧ

unsigned long


«

росток
«
ГСЧ

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;


/i> квантили нормального распределения


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;


/i> вид интегрируемой функции


0 —

постоянная


1 —

линейная


2 —

квадратичная

mcres.fun_type=Read1double(«fun_kind.txt»);


/i> вид системы ограничений


0

отсутствуют (весь параллелепипед)


1 —

линейные


2 —

квадратичное


3

сплайн
поверхности

mcres.con_type=Read1double(«con_type.txt»);


/i> загрузка парам е тров интегрируемой функции


switch

(mcres.fun_type)


{


case

fun_A=new
matd(«fun_A.txt»);


case

fun_b=new
matd(«fun_b.txt»);


case

fun_cd=Read1double(«fun_c.txt»);


}


/i> загрузка параметров ограничений


switch

(mcres.con_type)


{


case

сплайн — поверхности


/i> верхняя

xyz_top=
new

matd(«xyz_top.txt»);


/i> нижняя

xyz_bottom=
new

matd(«xyz_bottom.txt»);


/i> двумерная интерполяция

sp_top=
new

spl2d(xyz_top);

sp_bottom=
new

spl2d(xyz_bottom);


break;


case

квадратичная функция ограничений

con_U=
new

matd(«con_U.txt»);

con_v=
new

matd(«con_v.txt»);

con_wd=Read1double(«con_w.txt»);


break;


case

линейные ограничения

con_b=
new

matd(«con_b.txt»); con_A=new
matd(«con_A.txt»);


}


/i> объемлющий параллелепипед

a_int=
new

matd(«con_xmin.txt»);

b_int=
new

matd(«con_xmax.txt»);


/i> разность границ параллелепипеда

ba_int=
new

matd;

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


/i> аргумент интегрируемой функции

x_int=
new

matd(d_int,1);


/i> объем объемлющего параллелепипеда

mcres.V0_int=1;


for

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


{

/b> (_p(ba_int,j,1) <= 0)


{

Box(«Нижняя граница объемлющего параллелепипеда выше верхней для \

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


goto

clean_exit;


}

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


}


/i> начальный объем выборки

mcres.n1_int=10000;


/i> основной цикл для достижения заданной точности


/i> число итераций, потребовавшихся для достижения заданной то ч ности

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();


/i> генерация случайных точек и накопление суммы

sum1_int=0; sum2_int=0;

mcres.in_G_int=0;

PSChunk=
long

(mcres.n1_int/50.0);


/i> запуск ГСЧ

r=mcres.rng_seed;


for

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


{

b> int (r/m_rng);

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

/b> (r > d_rng) r=r-d_rng;


}


for

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


{


/i> случайный вектор


for

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


{


/i> случайное число

b> int (r/m_rng);

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

/b> (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;


}


/i> прогресс

/b> (!(i % PSChunk))


{

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

WaitForm->Refresh();


}


/i> проверка ограничения

con_ok=1;


switch

(mcres.con_type)

{


case

сплайн — поверхности

/b> ((_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

квадратичная функция ограничений

con_sum=0;


for

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


for

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

/b> (_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; <= d_int; ii++)

/b> (_p(con_v,ii,1) != 0)

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

/b> (con_sum > con_wd) con_ok=0;


break ;


case

линейная функция ограничений


for

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


{

con_sum=0;


for

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

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

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


}


}

fu_int=0;

/b> (con_ok != 0)


{

mcres.in_G_int++;


/i> точки графика

/b> (d_int==3)

/b> (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);


}


/i> значение интегрируемой функции


switch

(mcres.fun_type)


{


case

квадратичный член


for

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


for

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

/b> (_p(fun_A,ii,jj) != 0)

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

case
линейный член


for

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

/b> (_p(fun_b,ii,1) != 0)

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

case
постоянная

fu_int += fun_cd;


}


}

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


}


/i> оценка мат. ожидания и дисперсии

mcres.f1_int=sum1_int/mcres.n1_int;

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


/i> расчет погрешности

/b> (mcres.relok==0)


{


/i> абсолютная погрешность

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


}


else


{


/i> относительная погрешность

/b> (mcres.f1_int!=0)


{

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


}


else


{


/i> форма результатов

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();

/b> (mcres.deltar < mcres.dlt_int)


{


/i> точность достаточна

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

;


}


/i> вычисление нового объема выборки

/b> (mcres.relok==0)


{


/i> абс . погрешность

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


}


else


{


/i> отн . погрешность

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


}


/i> корректировка объема выборки в большую сторону


/i> для сокращения числа итер а ций

mcres.n1_int=1.2*mcres.n1_int;


/i> минимальный объем выборки

/b> (mcres.n1_int < 1000) mcres.n1_int=1000;

}
коне ц основно г о цикла

WaitForm->Close();


/i> график

/b> (d_int==3)


{

/b> (mcres.in_G_int==0)


{


/i> множество точек пусто

Zero_File(«xyz.txt»);


}


else

/b> (mcres.in_G_int < xyz->nl)


{


/i> точек не набралось, чтобы заполнить матрицу

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


{


/i> вся матрица заполнена

xyz->txprint(«xyz.txt»);


}


}

конец d_int==3

clean_exit:


/i> очистка памяти

/b> (d_int==3) delete xyz;


switch

(mcres.fun_type)


{


case

delete
fun_A;


case

delete
fun_b;


}


switch

(mcres.con_type)


{


case

delete
xyz_top,xyz_bottom,sp_top,sp_bottom; break
;


case

delete
con_U,con_v; break
;


case

delete
con_b,con_A;


}


delete

a_int,b_int,ba_int,x_int;


}

//integral

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


struct

mcres_struct


{


/i> индикатор относительной погрешности


int

relok;


/i> целевая погрешность


double

dlt_int;


/i> достигнутая погрешность


double

deltar;


/i> доверительная вероятность


double

omega_int;


/i> квантиль норм. распределения


double

z_int;


«

росток
«
ГСЧ


unsigned

long
rng_seed;


?

число итераций, потребовавшихся для достижения заданной точн
о
сти


unsigned

n_ite;


/i> объем выборки на последней итерации


unsigned

long
n1_int;


/i> число точек попавших в область интегрирования


unsigned

in_G_int;


/i> интеграл


double

inte_int;


/i> объем объемлющего параллелепипеда


double

V0_int;


/i> выборочное среднее


double

f1_int;


/i> выборочная дисперсия


double

vari_int;


/i> время начала счета

time_t t_start;


/i> время окончания счета

time_t t_end;


/i> продолжительность вычисления интеграла

time_t t_calc;


/i> вид интегрируемой функции


int

fun_type;


/i> вид системы огрничений


int

con_type;


}

; mcres_struct

ПРИЛОЖЕНИЕ 2

ТЕСТОВЫЕ ПРИМЕРЫ

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

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

Приближенное значение найдено для целевой абсолютной погрешности 0.00001.

Погрешность: 0.000034416630896 или 0.014749984670 %.

Примеры 2-10 Объемы многомерных шаров

Точные и приближенные объемы многомерных шаров приведены в следующей таблице.

Объем

точный Источник [5], с. 287.

Объем приближенный Вычислено в Maple (20 значащих цифр).

Оценка CarloS Расчет с целевой относительной погрешностью 2%

Относительная погрешность, %

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

2.5501640398773454440

2.55096385956571

.31363460384e-1