1 // Copyright (c) 1997-1999 Matra Datavision
2 // Copyright (c) 1999-2014 OPEN CASCADE SAS
4 // This file is part of Open CASCADE Technology software library.
6 // This library is free software; you can redistribute it and/or modify it under
7 // the terms of the GNU Lesser General Public License version 2.1 as published
8 // by the Free Software Foundation, with special exception defined in the file
9 // OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT
10 // distribution for complete text of the license and disclaimer of any warranty.
12 // Alternatively, this file may be used under the terms of Open CASCADE
13 // commercial license or contractual agreement.
16 Ce programme permet le calcul de l'integrale nieme d'une fonction. La
17 resolution peut etre vue comme l'integrale de l'integrale de ....(n fois), ce
18 qui necessite l'emploi de fonction recursive.
20 Une integrale simple de Gauss est la somme de la fonction donnee sur les points
21 de Gauss affectee du poids de Gauss correspondant.
23 I = somme(F(GaussPt) * GaussW) sur i = 1 jusqu'a l'ordre d'integration.
25 Une integrale multiple est une somme recursive sur toutes les variables.
27 Tout le calcul est effectue dans la methode recursive_iteration appelee n fois.
28 (n = nombre de variables)
30 Cette methode fait une sommation sur toutes les variables, sur tous les ordres
31 d'integration correspondants de la valeur de la fonction aux points et poids
37 #define No_Standard_RangeError
38 #define No_Standard_OutOfRange
39 #define No_Standard_DimensionError
44 #include <math_GaussMultipleIntegration.hxx>
45 #include <math_IntegerVector.hxx>
46 #include <math_Matrix.hxx>
47 #include <math_MultipleVarFunction.hxx>
48 #include <math_Vector.hxx>
49 #include <StdFail_NotDone.hxx>
51 class IntegrationFunction {
53 // Cette classe sert a conserver dans les champs les valeurs de tests de
54 // la fonction recursive.(En particulier le pointeur sur la fonction F,
55 // heritant d'une classe abstraite, non traite en cdl.)
57 math_MultipleVarFunction *Fsav;
58 math_IntegerVector Ordsav;
59 Standard_Integer NVarsav;
62 math_Matrix GaussPoint;
63 math_Matrix GaussWeight;
65 Standard_Boolean Done;
68 IntegrationFunction(math_MultipleVarFunction& F,
69 const Standard_Integer maxsav, const Standard_Integer NVar,
70 const math_IntegerVector& Ord,
71 const math_Vector& Lowsav,
72 const math_Vector& Uppsav);
74 Standard_Real Value() ;
75 Standard_Boolean IsDone() const;
76 Standard_Boolean recursive_iteration(Standard_Integer& n, math_IntegerVector& inc);
79 IntegrationFunction::IntegrationFunction(math_MultipleVarFunction& F,
80 const Standard_Integer maxsav, const Standard_Integer NVar,
81 const math_IntegerVector& Ord,
82 const math_Vector& Lowsav,const math_Vector& Uppsav):
86 GaussPoint(1, NVar, 1, maxsav),
87 GaussWeight(1, NVar, 1, maxsav)
90 Standard_Integer i, k;
91 math_IntegerVector inc(1, NVar);
96 Done = Standard_False;
98 //Recuperation des points et poids de Gauss dans le fichier GaussPoints
99 for (i =1; i<= NVarsav; i++) {
100 xm(i) = 0.5*(Lowsav(i) + Uppsav(i));
101 xr(i) = 0.5*(Uppsav(i) - Lowsav(i));
102 math_Vector GP(1,Ordsav(i)), GW(1,Ordsav(i));
103 math::GaussPoints(Ordsav(i),GP);
104 math::GaussWeights(Ordsav(i),GW);
105 for (k = 1; k <= Ordsav(i); k++) {
106 GaussPoint(i,k) = GP(k); //kieme point et poids de
107 GaussWeight(i,k) = GW(k);//Gauss de la variable i.
111 Standard_Integer Iterdeb = 1;
112 Standard_Boolean recur = recursive_iteration(Iterdeb, inc);
114 //On ramene l'integration a la bonne echelle.
115 for ( i = 1; i <= NVarsav; i++) {
118 Done = Standard_True;
122 Standard_Real IntegrationFunction::Value() {
126 Standard_Boolean IntegrationFunction::IsDone() const{
130 Standard_Boolean IntegrationFunction::recursive_iteration(Standard_Integer& n,
131 math_IntegerVector& inc) {
133 // Termination criterium :
134 // Calcul de la valeur de la fonction aux points de Gauss fixes:
136 if (n == (NVarsav+1)) {
137 math_Vector dx(1, NVarsav);
139 for ( j = 1; j <= NVarsav; j++) {
140 dx(j) = xr(j)* GaussPoint(j, inc(j));
143 Standard_Boolean Ok = Fsav->Value(xm + dx, F1);
144 if (!Ok) {return Standard_False;};
145 Standard_Real Interm = 1;
146 for ( j = 1; j <= NVarsav; j++) {
147 Interm *= GaussWeight(j, inc(j));
150 return Standard_True;
153 // Calcul recursif, pour chaque variable n de la valeur de la fonction aux
154 // Ordsav(n) points de Gauss.
156 Standard_Boolean OK=Standard_False;
157 for (inc(n) = 1; inc(n) <= Ordsav(n); inc(n)++) {
159 OK = recursive_iteration(local, inc);
164 math_GaussMultipleIntegration::
165 math_GaussMultipleIntegration(math_MultipleVarFunction& F,
166 const math_Vector& Lower,
167 const math_Vector& Upper,
168 const math_IntegerVector& Order)
170 Standard_Integer MaxOrder = math::GaussPointsMax();
172 Standard_Integer i, max =0;
173 Standard_Integer NVar = F.NbVariables();
174 math_IntegerVector Ord(1, NVar);
175 math_Vector Lowsav(1, NVar);
176 math_Vector Uppsav(1, NVar);
180 Done = Standard_False;
181 for (i = 1; i <= NVar; i++) {
182 if (Order(i) > MaxOrder) {
188 if(Ord(i) >= max) {max = Ord(i);}
191 //Calcul de l'integrale aux points de Gauss.
192 // Il s agit d une somme multiple sur le domaine [Lower, Upper].
194 IntegrationFunction Func(F, max, NVar, Ord, Lowsav, Uppsav);
197 Done = Standard_True;
201 void math_GaussMultipleIntegration::Dump(Standard_OStream& o) const {
203 o <<"math_GaussMultipleIntegration ";
205 o << " Status = Done \n";
206 o << " Integration value = " << Val <<"\n";
209 o << "Status = not Done \n";