Листинг 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] Источник [5], с. 287. [2] Вычислено в Maple (20 значащих цифр). [3] Расчет с целевой относительной погрешностью 2%
Популярное: Как построить свою речь (словесное оформление):
При подготовке публичного выступления перед оратором возникает вопрос, как лучше словесно оформить свою... Как вы ведете себя при стрессе?: Вы можете самостоятельно управлять стрессом! Каждый из нас имеет право и возможность уменьшить его воздействие на нас... Почему двоичная система счисления так распространена?: Каждая цифра должна быть как-то представлена на физическом носителе... ©2015-2024 megaobuchalka.ru Все материалы представленные на сайте исключительно с целью ознакомления читателями и не преследуют коммерческих целей или нарушение авторских прав. (166)
|
Почему 1285321 студент выбрали МегаОбучалку... Система поиска информации Мобильная версия сайта Удобная навигация Нет шокирующей рекламы |