Математика: Вычисление интегралов методом Монте-Карло, Курсовая работа

ПРИЛОЖЕНИЕ 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;

  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