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
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.
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.
13 I = somme(F(GaussPt) * GaussW) sur i = 1 jusqu'a l'ordre d'integration.
15 Une integrale multiple est une somme recursive sur toutes les variables.
17 Tout le calcul est effectue dans la methode recursive_iteration appelee n fois.
18 (n = nombre de variables)
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
27 #define No_Standard_RangeError
28 #define No_Standard_OutOfRange
29 #define No_Standard_DimensionError
32 #include <math_GaussMultipleIntegration.ixx>
34 #include <math_Matrix.hxx>
35 #include <math_Vector.hxx>
36 #include <math_IntegerVector.hxx>
37 #include <math_MultipleVarFunction.hxx>
39 class IntegrationFunction {
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.)
45 math_MultipleVarFunction *Fsav;
46 math_IntegerVector Ordsav;
47 Standard_Integer NVarsav;
50 math_Matrix GaussPoint;
51 math_Matrix GaussWeight;
53 Standard_Boolean Done;
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);
62 Standard_Real Value() ;
63 Standard_Boolean IsDone() const;
64 Standard_Boolean recursive_iteration(Standard_Integer& n, math_IntegerVector& inc);
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):
74 GaussPoint(1, NVar, 1, maxsav),
75 GaussWeight(1, NVar, 1, maxsav)
78 Standard_Integer i, k;
79 math_IntegerVector inc(1, NVar);
84 Done = Standard_False;
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.
99 Standard_Integer Iterdeb = 1;
100 Standard_Boolean recur = recursive_iteration(Iterdeb, inc);
102 //On ramene l'integration a la bonne echelle.
103 for ( i = 1; i <= NVarsav; i++) {
106 Done = Standard_True;
110 Standard_Real IntegrationFunction::Value() {
114 Standard_Boolean IntegrationFunction::IsDone() const{
118 Standard_Boolean IntegrationFunction::recursive_iteration(Standard_Integer& n,
119 math_IntegerVector& inc) {
121 // Termination criterium :
122 // Calcul de la valeur de la fonction aux points de Gauss fixes:
124 if (n == (NVarsav+1)) {
125 math_Vector dx(1, NVarsav);
127 for ( j = 1; j <= NVarsav; j++) {
128 dx(j) = xr(j)* GaussPoint(j, inc(j));
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));
138 return Standard_True;
141 // Calcul recursif, pour chaque variable n de la valeur de la fonction aux
142 // Ordsav(n) points de Gauss.
144 Standard_Boolean OK=Standard_False;
145 for (inc(n) = 1; inc(n) <= Ordsav(n); inc(n)++) {
147 OK = recursive_iteration(local, inc);
152 math_GaussMultipleIntegration::
153 math_GaussMultipleIntegration(math_MultipleVarFunction& F,
154 const math_Vector& Lower,
155 const math_Vector& Upper,
156 const math_IntegerVector& Order)
158 Standard_Integer MaxOrder = math::GaussPointsMax();
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);
168 Done = Standard_False;
169 for (i = 1; i <= NVar; i++) {
170 if (Order(i) > MaxOrder) {
176 if(Ord(i) >= max) {max = Ord(i);}
179 //Calcul de l'integrale aux points de Gauss.
180 // Il s agit d une somme multiple sur le domaine [Lower, Upper].
182 IntegrationFunction Func(F, max, NVar, Ord, Lowsav, Uppsav);
185 Done = Standard_True;
189 void math_GaussMultipleIntegration::Dump(Standard_OStream& o) const {
191 o <<"math_GaussMultipleIntegration ";
193 o << " Status = Done \n";
194 o << " Integration value = " << Val <<"\n";
197 o << "Status = not Done \n";