МАТЕМАТИЧЕСКОЕ ОБОСНОВАНИЕ АЛГОРИТМА ВЫЧИСЛЕНИЯ ИНТЕГРАЛА

Поможем в ✍️ написании учебной работы

Имя

Поможем с курсовой, контрольной, дипломной, рефератом, отчетом по практике, научно-исследовательской и любой другой работой

Нажимая кнопку "Продолжить", я принимаю политику конфиденциальности

Задание 12 Вычисление интегралов методом Монте – Карло.

Цель:

1) Реализация генератора случайных чисел для метода Монте – Карло.

2) Сравнение равномерного распределения и специально разработанного.

3) Вычисление тестового многомерного интеграла в сложной области.

Продукт:

1) Программный код в виде функции на языке С++ или Fortran .

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

3) Обзор использованной литературы.

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

МАТЕМАТИЧЕСКОЕ ОБОСНОВАНИЕ АЛГОРИТМА ВЫЧИСЛЕНИЯ ИНТЕГРАЛА

Принцип работы метода Монте – Карло

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

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

Сущность метода состоит в том, что в задачу вводят случайную величину , изменяющуюся по какому то правилу . Случайную величину выбирают таким образом, чтобы искомая в задаче величина стала математическим ожидание от , то есть .

Таким образом, искомая величина определяется лишь теоретически. Чтобы найти ее численно необходимо воспользоваться статистическими методами. То есть необходимо взять выборку случайных чисел объемом . Затем необходимо вычислить выборочное среднее варианта случайной величины по формуле:

. (1)

Вычисленное выборочное среднее принимают за приближенное значение .

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

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

Сплайн – интерполяция.

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

Под сплайном (от англ. spline - планка, рейка) обычно понимают агрегатную функцию, совпадающую с функциями более простой природы на каждом элементе разбиения своей области определения. Сплайн – функция имеет следующий вид:

. (15)

Исходные данные представляют собой троек точек .

Коэффициенты и определяются из системы:

, (16)

где ,

Алгоритм расчета интеграла

Реализованный алгоритм включает следующие шаги:

1) выбирается начальное значение , разыгрываются случайные векторы из и определяются и ;

2) в зависимости от вида погрешности (абсолютная, относительная) определяется достигнутая погрешность; если она меньше целевой, вычисление прерывается;

3) по формулам (13) или (14) вычисляется новый объем выборки;

4) объем выборки увеличивается на 20%

5) переход к шагу 1;

6) конец.

ЗАКЛЮЧЕНИЕ

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

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

Листинг 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 Объемы многомерных шаров

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

	Объем точный[1]	Объем приближенный[2]	Оценка 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	.31363460384e-1