src/math/math_GaussSingleIntegration.cxx

   1 // File math_GaussSingleIntegration.cxx
   2
   3 /*
   4 Par Gauss le calcul d'une integrale simple se transforme en sommation des
   5 valeurs de la fonction donnee aux <Order> points de Gauss affectee des poids
   6 de Gauss.
   7
   8 Les points et poids de Gauss sont stockes dans GaussPoints.cxx.
   9 Les points sont compris entre les valeurs -1 et +1, ce qui necessite un
  10 changement de variable pour les faire varier dans l'intervalle [Lower, Upper].
  11
  12
  13 On veut calculer Integrale( f(u)* du) entre a et b.
  14
  15
  16 Etapes du calcul:
  17
  18 1- calcul de la fonction au ieme point de Gauss (apres changement de variable).
  19
  20 2- multiplication de cette valeur par le ieme poids de Gauss.
  21
  22 3- sommation de toutes ces valeurs.
  23
  24 4- retour a l'intervalle [Lower, Upper] de notre integrale.
  25
  26 */
  27
  28 //#ifndef DEB
  29 #define No_Standard_RangeError
  30 #define No_Standard_OutOfRange
  31 #define No_Standard_DimensionError
  32 //#endif
  33
  34 #include <math_GaussSingleIntegration.ixx>
  35
  36 #include <math.hxx>
  37 #include <math_Vector.hxx>
  38 #include <math_Function.hxx>
  39
  40 math_GaussSingleIntegration::math_GaussSingleIntegration() : Done(Standard_False)
  41 {
  42 }
  43
  44 math_GaussSingleIntegration::
  45                  math_GaussSingleIntegration(math_Function& F,
  46                                              const Standard_Real Lower,
  47                                              const Standard_Real Upper,
  48                                              const Standard_Integer Order)
  49 {
  50   Standard_Integer theOrder = Min(math::GaussPointsMax(), Order);
  51   Perform(F, Lower, Upper, theOrder);
  52 }
  53
  54 math_GaussSingleIntegration::
  55                  math_GaussSingleIntegration(math_Function& F,
  56                                              const Standard_Real Lower,
  57                                              const Standard_Real Upper,
  58                                              const Standard_Integer Order,
  59                                              const Standard_Real Tol)
  60 {
  61   Standard_Integer theOrder = Min(math::GaussPointsMax(), Order);
  62
  63   const Standard_Integer IterMax = 13;   // Max number of iteration
  64   Standard_Integer NIter = 1;            // current number of iteration
  65   Standard_Integer NbInterval = 1;       // current number of subintervals
  66   Standard_Real    dU,OldLen,Len;
  67
  68   Perform(F, Lower, Upper, theOrder);
  69   Len = Val;
  70   do {
  71     OldLen = Len;
  72     Len = 0.;
  73     NbInterval *= 2;
  74     dU = (Upper-Lower)/NbInterval;
  75     for (Standard_Integer i=1; i<=NbInterval; i++) {
  76       Perform(F, Lower+(i-1)*dU, Lower+i*dU, theOrder);
  77       if (!Done) return;
  78       Len += Val;
  79     }
  80     NIter++;
  81   }
  82   while (fabs(OldLen-Len) > Tol && NIter <= IterMax);
  83
  84   Val = Len;
  85 }
  86
  87 void math_GaussSingleIntegration::Perform(math_Function& F,
  88                                           const Standard_Real Lower,
  89                                           const Standard_Real Upper,
  90                                           const Standard_Integer Order)
  91 {
  92   Standard_Real xr, xm, dx;
  93   Standard_Integer j;
  94   Standard_Real F1, F2;
  95   Standard_Boolean Ok1;
  96   math_Vector GaussP(1, Order);
  97   math_Vector GaussW(1, Order);
  98   Done = Standard_False;
  99
 100 //Recuperation des points de Gauss dans le fichier GaussPoints.
 101   math::GaussPoints(Order,GaussP);
 102   math::GaussWeights(Order,GaussW);
 103
 104 // Calcul de l'integrale aux points de Gauss :
 105
 106 // Changement de variable pour la mise a l'echelle [Lower, Upper] :
 107   xm = 0.5*(Upper + Lower);
 108   xr = 0.5*(Upper - Lower);
 109   Val = 0.;
 110
 111   Standard_Integer ind = Order/2, ind1 = (Order+1)/2;
 112   if(ind1 > ind) { // odder case
 113     Ok1 = F.Value(xm, Val);
 114     if (!Ok1) return;
 115     Val *= GaussW(ind1);
 116   }
 117 // Sommation sur tous les points de Gauss: avec utilisation de la symetrie.
 118   for (j = 1; j <= ind; j++) {
 119     dx = xr*GaussP(j);
 120     Ok1 = F.Value(xm-dx, F1);
 121     if(!Ok1) return;
 122     Ok1 = F.Value(xm+dx, F2);
 123     if(!Ok1) return;
 124     // Multiplication par les poids de Gauss.
 125     Standard_Real FT = F1+F2;
 126     Val += GaussW(j)*FT;
 127   }
 128   // Mise a l'echelle de l'intervalle [Lower, Upper]
 129   Val *= xr;
 130   Done = Standard_True;
 131 }
 132
 133 void math_GaussSingleIntegration::Dump(Standard_OStream& o) const {
 134
 135   o <<"math_GaussSingleIntegration ";
 136    if (Done) {
 137      o << " Status = Done \n";
 138      o << "Integration Value = " << Val<<"\n";
 139    }
 140    else {
 141      o << "Status = not Done \n";
 142    }
 143 }