6-мерный интеграл по трапеции на Фортране с использованием Фортрана 90

Вы можете посмотреть главу об интеграции Монте-Карло в Научной библиотеке GNU (GSL). Это и библиотека, и, поскольку это открытый исходный код, исходный код, который вы можете изучить.

Посмотрите раздел 4.6 числовых рецептов для C.

1. Первый шаг — уменьшить проблему с помощью симметрии и аналитических зависимостей.

2. Шаг второй — связать решение следующим образом:

   f2(x2,x3,..,x6) = Integrate(f(x,x2,x3..,x6),x,1,x1end)
   f3(x3,x4,..,x6) = Integrate(f2(x,x3,..,x6),x,1,x2end)
   f4(x4,..,x6) = ...
   
   f6(x6) = Integrate(I4(x,x6),x,1,x5end)        
   result = Integrate(f6(x),x,1,x6end)
   

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

Процедура интеграции, которую я использую, называется «QuadMo», написанной Люком Мо примерно в 1970 году. Я сделал ее рекурсивной и поместил в модуль. QuadMo уточняет сетку, необходимую для получения требуемой точности интегрирования. Вот программа, которая вычисляет n-мерный интеграл, используя QuadMo.

Вот проверка программы с использованием Гаусса с центром в 0,5 и SD 0,1 во всех измерениях для nDim до 6 с использованием компиляции G95. Запускается за пару секунд.

  nDim ans       expected       nlvl
     1 0.249     0.251      2
     2 6.185E-02 6.283E-02  2  2
     3 1.538E-02 1.575E-02  2  2  2
     4 3.826E-03 3.948E-03  2  2  2  2
     5 9.514E-04 9.896E-04  2  2  2  2  2
     6 2.366E-04 2.481E-04  2  2  2  2  2  2

Вот код:

!=======================================================================
      module QuadMo_MOD
      implicit none
      integer::QuadMo_MinLvl=6,QuadMo_MaxLvl=24
      integer,dimension(:),allocatable::QuadMo_nlvlk
      real*8::QuadMo_Tol=1d-5
      real*8,save,dimension(:),allocatable::thet
      integer,save::nDim

      abstract interface
        function QuadMoFunct_interface(thet,k)
           real*8::QuadMoFunct_interface
           real*8,intent(in)::thet
           integer,intent(in),optional::k
        end function
      end interface

      abstract interface
        function MultIntFunc_interface(thet)
           real*8::MultIntFunc_interface
           real*8,dimension(:),intent(in)::thet
        end function
      end interface

      procedure(MultIntFunc_interface),pointer :: stored_func => null()

      contains

!----------------------------------------------------------------------
      recursive function quadMoMult(funct,lower,upper,k) result(ans) 
  ! very powerful integration routine written by Luke Mo
  !  then at the Stanford Linear Accelerator Center circa 1970
  ! QuadMo_Eps is error tolerance  
  ! QuadMo_MinLvl determines initial grid of 2**(MinLvl+1) + 1 points
  !  to avoid missing a narrow peak, this may need to be increased.  
  ! QuadMo_Nlvl returns number of subinterval refinements required beyond
  ! QuadMo_MaxLvl

  ! Modified by making recursive and adding argument k
  !  for multiple integrals ([email protected])
      real*8::ans
      procedure(QuadMoFunct_interface) :: funct
      real*8,intent(in)::lower,upper
      integer,intent(in),optional::k
      real*8::Middle,Left,Right,eps,est,fLeft,fMiddle,fRight
     & ,fml,fmr,rombrg,coef,estl,estr,estint,area,abarea
      real*8::valint(50,2), Middlex(50), Rightx(50), fmx(50), frx(50)
     & ,fmrx(50), estrx(50), epsx(50)
      integer retrn(50),i,level
      level = 0
      QuadMo_nlvlk(k) = 0
      abarea = 0
      Left = lower
      Right = upper
      if(present(k))then
         fLeft = funct(Left,k)
         fMiddle = funct((Left+Right)/2,k)
         fRight = funct(Right,k)
      else
         fLeft = funct(Left)
         fMiddle = funct((Left+Right)/2)
         fRight = funct(Right)
      endif
      est = 0
      eps = QuadMo_Tol
  100 level = level+1
      Middle = (Left+Right)/2
      coef = Right-Left
      if(coef.ne.0) go to 150
      rombrg = est
      go to 300
  150 continue
      if(present(k))then
         fml = funct((Left+Middle)/2,k)
         fmr = funct((Middle+Right)/2,k)
      else
         fml = funct((Left+Middle)/2)
         fmr = funct((Middle+Right)/2)
      endif
      estl = (fLeft+4*fml+fMiddle)*coef
      estr = (fMiddle+4*fmr+fRight)*coef
      estint = estl+estr
      area= abs(estl)+ abs(estr)
      abarea=area+abarea- abs(est)
      if(level.ne.QuadMo_MaxLvl) go to 200
      QuadMo_nlvlk(k) = QuadMo_nlvlk(k)+1
      rombrg = estint
      go to 300
  200 if(( abs(est-estint).gt.(eps*abarea)).or.
     1(level.lt.QuadMo_MinLvl))  go to 400
      rombrg = (16*estint-est)/15
  300 level = level-1
      i = retrn(level)
      valint(level, i) = rombrg
      go to (500, 600), i
  400 retrn(level) = 1
      Middlex(level) = Middle
      Rightx(level) = Right
      fmx(level) = fMiddle
      fmrx(level) = fmr
      frx(level) = fRight
      estrx(level) = estr
      epsx(level) = eps
      eps = eps/1.4d0
      Right = Middle
      fRight = fMiddle
      fMiddle = fml
      est = estl
      go to 100
  500 retrn(level) = 2
      Left = Middlex(level)
      Right = Rightx(level)
      fLeft = fmx(level)
      fMiddle = fmrx(level)
      fRight = frx(level)
      est = estrx(level)
      eps = epsx(level)
      go to 100
  600 rombrg = valint(level,1)+valint(level,2)
      if(level.gt.1) go to 300
      ans  = rombrg /12
      end function quadMoMult
!-----------------------------------------------------------------------
      recursive function MultInt(k,func) result(ans)
      ! MultInt(nDim,func) returns multi-dimensional integral from 0 to 1
      !  in all dimensions of function func
      !  variable QuadMo_Mod: nDim needs to be set initially to number of dimensions      
      procedure(MultIntFunc_interface) :: func
      real*8::ans
      integer,intent(in)::k

      stored_func => func

      if(k.eq.nDim)then
         if(allocated(thet))deallocate(thet)
         allocate (thet(nDim))
         if(allocated(QuadMo_nlvlk))deallocate(QuadMo_nlvlk)
         allocate(QuadMo_nlvlk(nDim))
      endif

      if(k.eq.0)then
         ans=func(thet)
         return
      else
         ans=QuadMoMult(MultIntegrand,0d0,1d0,k)
      endif
      end function MultInt
!-----------------------------------------------------------------------
      recursive function MultIntegrand(thetARG,k) result(ans)
      real*8::ans
      real*8,intent(in)::thetARG
      integer,optional,intent(in)::k

      if(present(k))then
         thet(k)=thetARG
      else
         write(*,*)'MultIntegrand: not expected, k not present!'
         stop
      endif
      ans=MultInt(k-1,stored_func)
      end function MultIntegrand
!-----------------------------------------------------------------------
      end module QuadMo_MOD
!=======================================================================
      module test_MOD
      use QuadMo_MOD
      implicit none

      contains

!----------------------------------------------------------------------- 
      real*8 function func(thet) ! multidimensional function
      ! this is the function defined in nDim dimensions
      !  in this case a Gaussian centered at 0.5 with SD 0.1
      real*8,intent(in),dimension(:)::thet
      func=exp(-sum(((thet-5d-1)/1d-1)
     & *((thet-5d-1)/1d-1))/2)
      end function func
!----------------------------------------------------------------------- 
      end module test_MOD 
!=======================================================================
      ! test program to evaluate multiple integrals 
      use test_MOD
      implicit none 
      real*8::ans 

      ! these values are set for speed, not accuracy
      QuadMo_MinLvl=2
      QuadMo_MaxLvl=3
      QuadMo_Tol=1d-1

      write(*,*)'     nDim ans       expected       nlvl'
      do nDim=1,6
         ! expected answer is (0.1 sqrt(2pi))**nDim
         ans=MultInt(nDim,func)
         write(*,'(i10,2(1pg10.3),999(i3))')nDim,ans,(0.250663)**nDim
     &    ,QuadMo_nlvlk
      enddo
      end
!-----------------------------------------------------------------------

double MultInt(int k);
double MultIntegrand(double thetARG, int k);
double quadMoMult(double(*funct)(double, int), double lower, double upper, int k);
double funkn(double *thet);

int QuadMo_MinLvl = 2;
int  QuadMo_MaxLvl = 3;
double QuadMo_Tol = 0.1;
int *QuadMo_nlvlk;
double *thet;
int nDim;


//double MultInt(int k, double(*func)(double *))
double MultInt(int k)
{

//MultInt(nDim, func) returns multi - dimensional integral from 0 to 1
//in all dimensions of function func
    double ans;

    if (k == 0) 
    {
        ans = funkn(thet);
    }
    else
    {
        ans = quadMoMult(MultIntegrand, 0.0, 1.0, k); //limits hardcoded here
    }
    return ans;
}


double MultIntegrand(double thetARG, int k)
{
    double ans;

    if (k > 0)
        thet[k] = thetARG;
    else
        printf("\n***MultIntegrand: not expected, k not present!***\n");

        //Recursive call
    //ans = MultInt(k - 1, func);
    ans = MultInt(k - 1);
    return ans;
}


double quadMoMult(double(*funct)(double, int), double lower, double upper, int k)
{
    //Integration routine written by Luke Mo
    //Stanford Linear Accelerator Center circa 1970
    //QuadMo_Eps is error tolerance
    //QuadMo_MinLvl determines initial grid of 2 * *(MinLvl + 1) + 1 points
    //to avoid missing a narrow peak, this may need to be increased.
    //QuadMo_Nlvl returns number of subinterval refinements required beyond
    //QuadMo_MaxLvl

    //Modified by making recursive and adding argument k
    //for multiple integrals([email protected])

    double ans;
    double Middle, Left, Right, eps, est, fLeft, fMiddle, fRight;
    double fml, fmr, rombrg, coef, estl, estr, estint, area, abarea;
    double valint[51][3], Middlex[51], Rightx[51], fmx[51], frx[51];  //Jack up arrays
    double fmrx[51], estrx[51], epsx[51];
    int retrn[51];
    int i, level;

    level = 0;
    QuadMo_nlvlk[k] = 0;
    abarea = 0.0;
    Left = lower;
    Right = upper;
    if (k > 0)
    {
        fLeft = funct(Left, k);
        fMiddle = funct((Left + Right) / 2, k);
        fRight = funct(Right, k);
    }
    else
    {
        fLeft = funct(Left,0);
        fMiddle = funct((Left + Right) / 2,0);
        fRight = funct(Right,0);
    }
    est = 0.0;
    eps = QuadMo_Tol;
l100:
    level = level + 1;
    Middle = (Left + Right) / 2;
    coef = Right - Left;
    if (coef != 0.0)
        goto l150;
    rombrg = est;
    goto l300;
l150:
    if (k > 0)
    {
        fml = funct((Left + Middle) / 2.0, k);
        fmr = funct((Middle + Right) / 2.0, k);
    }
    else
    {
        fml = funct((Left + Middle) / 2.0, 0);
        fmr = funct((Middle + Right) / 2.0, 0);
    }
    estl = (fLeft + 4 * fml + fMiddle)*coef;
    estr = (fMiddle + 4 * fmr + fRight)*coef;
    estint = estl + estr;
    area = abs(estl) + abs(estr);
    abarea = area + abarea - abs(est);
    if (level != QuadMo_MaxLvl)
        goto l200;
    QuadMo_nlvlk[k] = QuadMo_nlvlk[k] + 1;
    rombrg = estint;
    goto l300;
l200:
    if ((abs(est - estint) > (eps*abarea)) || (level < QuadMo_MinLvl))
        goto l400;
    rombrg = (16 * estint - est) / 15;
l300:
    level = level - 1;
    i = retrn[level];
    valint[level][i] = rombrg;
    if (i == 1)
        goto l500;
    if (i == 2)
        goto l600;
l400:
    retrn[level] = 1;
    Middlex[level] = Middle;
    Rightx[level] = Right;
    fmx[level] = fMiddle;
    fmrx[level] = fmr;
    frx[level] = fRight;
    estrx[level] = estr;
    epsx[level] = eps;
    eps = eps / 1.4;
    Right = Middle;
    fRight = fMiddle;
    fMiddle = fml;
    est = estl;
    goto l100;
l500:
    retrn[level] = 2;
    Left = Middlex[level];
    Right = Rightx[level];
    fLeft = fmx[level];
    fMiddle = fmrx[level];
    fRight = frx[level];
    est = estrx[level];
    eps = epsx[level];
    goto l100;
l600:
    rombrg = valint[level][1] + valint[level][2];
    if (level > 1)
        goto l300;
    ans = rombrg / 12.0;
    return ans;
}


double funkn(double *thet)
{
    //in this case a Gaussian centered at 0.5 with SD 0.1
    double *sm;
    double sum;
    sm = new double[nDim];
    sum = 0.0;

    for (int i = 1; i <= nDim; i++) 
    {
        sm[i] = (thet[i] - 0.5) / 0.1;
        sm[i] *= sm[i];
        sum = sum + sm[i];
    }
    return exp(-sum / 2.0);
}


int main() {

    double ans;

    printf("\nnDim ans       expected       nlvl\n");
    for (nDim = 1; nDim <= 6; nDim++)
    {
        //expected answer is(0.1 sqrt(2pi))**nDim
        QuadMo_nlvlk = new int[nDim + 1];  //array for x values
        thet = new double[nDim + 1];  //array for x values
        ans = MultInt(nDim);
        printf("\n %d %f %f ", nDim, ans, pow((0.250663),nDim));
        for (int j=1; j<=nDim;  j++)
            printf("   %d ", QuadMo_nlvlk[nDim]);
        printf("\n");
    }

    return 0;

}

Объявите соответствующие параметры глобально

int QuadMo_MinLvl = 2;
int  QuadMo_MaxLvl = 3;
double QuadMo_Tol = 0.1;
int *QuadMo_nlvlk;
double *thet;
int nDim;

Это кодирование намного яснее, чем запутанное устаревшее кодирование на Фортране, с некоторой настройкой интегральные пределы и допуски могут быть параметризованы!!

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