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