7fd59977 |
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 | |