Вычисление интегралов методом Монте-Карло (стр. 3 из 3)
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 |
[1] Источник [5], с. 287.
[2] Вычислено в Maple (20 значащих цифр).
[3] Расчет с целевой относительной погрешностью 2%