src/math/math_GaussMultipleIntegration.cxx

   1 //static const char* sccsid = "@(#)math_GaussMultipleIntegration.cxx    3.3 95/01/25"; // Do not delete this line. Used by sccs.
   2 // File math_GaussMultipleIntegration.cxx
   3
   4
   5 /*
   6 Ce programme permet le calcul de l'integrale nieme d'une fonction. La
   7 resolution peut etre vue comme l'integrale de l'integrale de ....(n fois), ce
   8 qui necessite l'emploi de fonction recursive.
   9
  10 Une integrale simple de Gauss est la somme de la fonction donnee sur les points
  11 de Gauss affectee du poids de Gauss correspondant.
  12
  13 I = somme(F(GaussPt) * GaussW)   sur i = 1 jusqu'a l'ordre d'integration.
  14
  15 Une integrale multiple est une somme recursive sur toutes les variables.
  16
  17 Tout le calcul est effectue dans la methode recursive_iteration appelee n fois.
  18  (n = nombre de variables)
  19
  20 Cette methode fait une sommation sur toutes les variables, sur tous les ordres
  21 d'integration correspondants de la valeur de la fonction aux points et poids
  22 de Gauss.
  23
  24 */
  25
  26 //#ifndef DEB
  27 #define No_Standard_RangeError
  28 #define No_Standard_OutOfRange
  29 #define No_Standard_DimensionError
  30 //#endif
  31
  32 #include <math_GaussMultipleIntegration.ixx>
  33 #include <math.hxx>
  34 #include <math_Matrix.hxx>
  35 #include <math_Vector.hxx>
  36 #include <math_IntegerVector.hxx>
  37 #include <math_MultipleVarFunction.hxx>
  38
  39 class IntegrationFunction {
  40
  41 // Cette classe sert a conserver dans les champs les valeurs de tests de
  42 // la fonction recursive.(En particulier le pointeur sur la fonction F,
  43 // heritant d'une classe abstraite, non traite en cdl.)
  44
  45   math_MultipleVarFunction *Fsav;
  46   math_IntegerVector Ordsav;
  47   Standard_Integer NVarsav;
  48   math_Vector xr;
  49   math_Vector xm;
  50   math_Matrix GaussPoint;
  51   math_Matrix GaussWeight;
  52   Standard_Real Val;
  53   Standard_Boolean Done;
  54
  55  public:
  56     IntegrationFunction(math_MultipleVarFunction& F,
  57                         const Standard_Integer maxsav, const Standard_Integer NVar,
  58                         const math_IntegerVector& Ord,
  59                         const math_Vector& Lowsav,
  60                         const math_Vector& Uppsav);
  61
  62     Standard_Real Value() ;
  63     Standard_Boolean IsDone() const;
  64     Standard_Boolean recursive_iteration(Standard_Integer& n, math_IntegerVector& inc);
  65 };
  66
  67  IntegrationFunction::IntegrationFunction(math_MultipleVarFunction& F,
  68                         const Standard_Integer maxsav, const Standard_Integer NVar,
  69                         const math_IntegerVector& Ord,
  70                         const math_Vector& Lowsav,const  math_Vector& Uppsav):
  71                         Ordsav(1, NVar),
  72                         xr(1, NVar),
  73                         xm(1, NVar),
  74                         GaussPoint(1, NVar, 1, maxsav),
  75                         GaussWeight(1, NVar, 1, maxsav)
  76 {
  77
  78 Standard_Integer i, k;
  79 math_IntegerVector inc(1, NVar);
  80 inc.Init(0);
  81 Fsav = &F;
  82 NVarsav = NVar;
  83 Ordsav = Ord;
  84 Done = Standard_False;
  85
  86 //Recuperation des points et poids de Gauss dans le fichier GaussPoints
  87   for (i =1; i<= NVarsav; i++) {
  88     xm(i) = 0.5*(Lowsav(i) + Uppsav(i));
  89     xr(i) = 0.5*(Uppsav(i) - Lowsav(i));
  90     math_Vector GP(1,Ordsav(i)), GW(1,Ordsav(i));
  91     math::GaussPoints(Ordsav(i),GP);
  92     math::GaussWeights(Ordsav(i),GW);
  93     for (k = 1; k <= Ordsav(i); k++) {
  94       GaussPoint(i,k)  = GP(k);  //kieme point et poids de
  95       GaussWeight(i,k) = GW(k);//Gauss de la variable i.
  96     }
  97   }
  98   Val = 0.0;
  99   Standard_Integer Iterdeb = 1;
 100   Standard_Boolean recur = recursive_iteration(Iterdeb, inc);
 101   if (recur) {
 102 //On ramene l'integration a la bonne echelle.
 103      for ( i = 1; i <= NVarsav; i++) {
 104        Val *= xr(i);
 105      }
 106      Done = Standard_True;
 107   }
 108 }
 109
 110 Standard_Real IntegrationFunction::Value() {
 111   return Val;
 112 }
 113
 114 Standard_Boolean IntegrationFunction::IsDone() const{
 115   return Done;
 116 }
 117
 118 Standard_Boolean IntegrationFunction::recursive_iteration(Standard_Integer& n,
 119                                    math_IntegerVector& inc) {
 120
 121 // Termination criterium :
 122 // Calcul de la valeur de la fonction aux points de Gauss fixes:
 123         int local;
 124    if (n == (NVarsav+1)) {
 125      math_Vector dx(1, NVarsav);
 126      Standard_Integer j ;
 127      for ( j = 1; j <= NVarsav; j++) {
 128        dx(j) = xr(j)* GaussPoint(j, inc(j));
 129      }
 130      Standard_Real F1;
 131      Standard_Boolean Ok = Fsav->Value(xm + dx, F1);
 132      if (!Ok) {return Standard_False;};
 133      Standard_Real Interm = 1;
 134      for ( j = 1; j <= NVarsav; j++) {
 135        Interm *= GaussWeight(j, inc(j));
 136      }
 137      Val += Interm*F1;
 138      return Standard_True;
 139    }
 140
 141    // Calcul recursif, pour chaque variable n de la valeur de la fonction aux
 142    // Ordsav(n) points de Gauss.
 143
 144    Standard_Boolean OK=Standard_False;
 145    for (inc(n) = 1; inc(n) <= Ordsav(n); inc(n)++) {
 146         local = n + 1;
 147      OK = recursive_iteration(local, inc);
 148    }
 149    return OK;
 150 }
 151
 152 math_GaussMultipleIntegration::
 153                 math_GaussMultipleIntegration(math_MultipleVarFunction& F,
 154                                               const math_Vector& Lower,
 155                                               const math_Vector& Upper,
 156                                               const math_IntegerVector& Order)
 157 {
 158   Standard_Integer MaxOrder = math::GaussPointsMax();
 159
 160   Standard_Integer i,  max =0;
 161   Standard_Integer NVar = F.NbVariables();
 162   math_IntegerVector Ord(1, NVar);
 163   math_Vector Lowsav(1, NVar);
 164   math_Vector Uppsav(1, NVar);
 165   Lowsav = Lower;
 166   Uppsav = Upper;
 167 //  Ord = Order;
 168   Done = Standard_False;
 169   for (i = 1; i <= NVar; i++) {
 170     if (Order(i) > MaxOrder) {
 171       Ord(i) = MaxOrder;
 172     }
 173     else {
 174       Ord(i) = Order(i);
 175     }
 176     if(Ord(i) >= max) {max = Ord(i);}
 177   }
 178
 179 //Calcul de l'integrale aux points de Gauss.
 180 // Il s agit d une somme multiple sur le domaine [Lower, Upper].
 181
 182   IntegrationFunction Func(F, max, NVar, Ord, Lowsav, Uppsav);
 183   if (Func.IsDone()) {
 184     Val = Func.Value();
 185     Done = Standard_True;
 186   }
 187 }
 188
 189 void math_GaussMultipleIntegration::Dump(Standard_OStream& o) const {
 190
 191   o <<"math_GaussMultipleIntegration ";
 192    if (Done) {
 193      o << " Status = Done \n";
 194      o << " Integration value = " << Val <<"\n";
 195    }
 196    else {
 197      o << "Status = not Done \n";
 198    }
 199 }
 200
 201
 202