1 // Created on: 2005-12-08
2 // Created by: Sergey KHROMOV
3 // Copyright (c) 2005-2014 OPEN CASCADE SAS
5 // This file is part of Open CASCADE Technology software library.
7 // This library is free software; you can redistribute it and/or modify it under
8 // the terms of the GNU Lesser General Public License version 2.1 as published
9 // by the Free Software Foundation, with special exception defined in the file
10 // OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT
11 // distribution for complete text of the license and disclaimer of any warranty.
13 // Alternatively, this file may be used under the terms of Open CASCADE
14 // commercial license or contractual agreement.
18 #include <math_Function.hxx>
19 #include <math_KronrodSingleIntegration.hxx>
20 #include <math_Vector.hxx>
21 #include <StdFail_NotDone.hxx>
22 #include <TColStd_SequenceOfReal.hxx>
24 //==========================================================================
25 //function : An empty constructor.
26 //==========================================================================
27 math_KronrodSingleIntegration::math_KronrodSingleIntegration() :
28 myIsDone(Standard_False),
36 //==========================================================================
37 //function : Constructor
39 //==========================================================================
41 math_KronrodSingleIntegration::math_KronrodSingleIntegration
42 ( math_Function &theFunction,
43 const Standard_Real theLower,
44 const Standard_Real theUpper,
45 const Standard_Integer theNbPnts):
46 myIsDone(Standard_False),
51 Perform(theFunction, theLower, theUpper, theNbPnts);
54 //==========================================================================
55 //function : Constructor
57 //==========================================================================
59 math_KronrodSingleIntegration::math_KronrodSingleIntegration
60 ( math_Function &theFunction,
61 const Standard_Real theLower,
62 const Standard_Real theUpper,
63 const Standard_Integer theNbPnts,
64 const Standard_Real theTolerance,
65 const Standard_Integer theMaxNbIter) :
66 myIsDone(Standard_False),
71 Perform(theFunction, theLower, theUpper, theNbPnts,
72 theTolerance, theMaxNbIter);
75 //==========================================================================
77 // Computation of the integral.
78 //==========================================================================
80 void math_KronrodSingleIntegration::Perform
81 ( math_Function &theFunction,
82 const Standard_Real theLower,
83 const Standard_Real theUpper,
84 const Standard_Integer theNbPnts)
86 //const Standard_Real aMinVol = Epsilon(1.);
87 const Standard_Real aPtol = 1.e-9;
91 myIsDone = Standard_False;
95 if(theUpper - theLower < aPtol) {
96 myIsDone = Standard_False;
100 // Get an odd value of number of initial points.
101 myNbPntsReached = (theNbPnts%2 == 0) ? theNbPnts + 1 : theNbPnts;
102 myErrorReached = RealLast();
104 Standard_Integer aNGauss = myNbPntsReached/2;
105 math_Vector aKronrodP(1, myNbPntsReached);
106 math_Vector aKronrodW(1, myNbPntsReached);
107 math_Vector aGaussP(1, aNGauss);
108 math_Vector aGaussW(1, aNGauss);
111 if (!math::KronrodPointsAndWeights(myNbPntsReached, aKronrodP, aKronrodW) ||
112 !math::OrderedGaussPointsAndWeights(aNGauss, aGaussP, aGaussW)) {
113 myIsDone = Standard_False;
117 myIsDone = GKRule(theFunction, theLower, theUpper, aGaussP, aGaussW, aKronrodP, aKronrodW,
118 myValue, myErrorReached);
120 if(!myIsDone) return;
122 //Standard_Real anAbsVal = Abs(myValue);
124 myAbsolutError = myErrorReached;
126 //if (anAbsVal > aMinVol)
127 //myErrorReached /= anAbsVal;
133 //=======================================================================
137 //=======================================================================
139 void math_KronrodSingleIntegration::Perform
140 ( math_Function &theFunction,
141 const Standard_Real theLower,
142 const Standard_Real theUpper,
143 const Standard_Integer theNbPnts,
144 const Standard_Real theTolerance,
145 const Standard_Integer theMaxNbIter)
147 Standard_Real aMinVol = Epsilon(1.);
150 // Check prerequisites.
152 theTolerance <= 0.) {
153 myIsDone = Standard_False;
156 // Get an odd value of number of initial points.
157 myNbPntsReached = (theNbPnts%2 == 0) ? theNbPnts + 1 : theNbPnts;
159 Standard_Integer aNGauss = myNbPntsReached/2;
160 math_Vector aKronrodP(1, myNbPntsReached);
161 math_Vector aKronrodW(1, myNbPntsReached);
162 math_Vector aGaussP(1, aNGauss);
163 math_Vector aGaussW(1, aNGauss);
165 if (!math::KronrodPointsAndWeights(myNbPntsReached, aKronrodP, aKronrodW) ||
166 !math::OrderedGaussPointsAndWeights(aNGauss, aGaussP, aGaussW)) {
167 myIsDone = Standard_False;
172 myIsDone = GKRule(theFunction, theLower, theUpper, aGaussP, aGaussW, aKronrodP, aKronrodW,
173 myValue, myErrorReached);
175 if(!myIsDone) return;
177 Standard_Real anAbsVal = Abs(myValue);
179 myAbsolutError = myErrorReached;
180 if (anAbsVal > aMinVol)
181 myErrorReached /= anAbsVal;
185 if(myErrorReached <= theTolerance) return;
186 if(myNbIterReached >= theMaxNbIter) return;
188 TColStd_SequenceOfReal anIntervals;
189 TColStd_SequenceOfReal anErrors;
190 TColStd_SequenceOfReal aValues;
192 anIntervals.Append(theLower);
193 anIntervals.Append(theUpper);
196 anErrors.Append(myAbsolutError);
197 aValues.Append(myValue);
199 Standard_Integer i, nint, nbints;
201 Standard_Real maxerr;
202 Standard_Integer count = 0;
204 while(myErrorReached > theTolerance && myNbIterReached < theMaxNbIter) {
205 //Searching interval with max error
206 nbints = anIntervals.Length() - 1;
209 for(i = 1; i <= nbints; ++i) {
210 if(anErrors(i) > maxerr) {
211 maxerr = anErrors(i);
216 Standard_Real a = anIntervals(nint);
217 Standard_Real b = anIntervals(nint+1);
218 Standard_Real c = 0.5*(a + b);
220 Standard_Real v1, v2, e1, e2;
222 myIsDone = GKRule(theFunction, a, c, aGaussP, aGaussW, aKronrodP, aKronrodW,
225 if(!myIsDone) return;
227 myIsDone = GKRule(theFunction, c, b, aGaussP, aGaussW, aKronrodP, aKronrodW,
230 if(!myIsDone) return;
234 Standard_Real deltav = v1 + v2 - aValues(nint);
236 if(Abs(deltav) <= Epsilon(Abs(myValue))) ++count;
238 Standard_Real deltae = e1 + e2 - anErrors(nint);
239 myAbsolutError += deltae;
240 if(myAbsolutError <= Epsilon(Abs(myValue))) ++count;
242 if(Abs(myValue) > aMinVol) myErrorReached = myAbsolutError/Abs(myValue);
243 else myErrorReached = myAbsolutError;
246 if(count > 50) return;
248 //Inserting new interval
250 anIntervals.InsertAfter(nint, c);
252 anErrors.InsertAfter(nint, e2);
254 aValues.InsertAfter(nint, v2);
260 //=======================================================================
263 //=======================================================================
265 Standard_Boolean math_KronrodSingleIntegration::GKRule(
266 math_Function &theFunction,
267 const Standard_Real theLower,
268 const Standard_Real theUpper,
269 const math_Vector& /*theGaussP*/,
270 const math_Vector& theGaussW,
271 const math_Vector& theKronrodP,
272 const math_Vector& theKronrodW,
273 Standard_Real& theValue,
274 Standard_Real& theError)
277 Standard_Boolean IsDone;
279 Standard_Integer aNKronrod = theKronrodP.Length();
281 Standard_Real aGaussVal;
282 Standard_Integer aNPnt2 = (aNKronrod + 1)/2;
288 math_Vector f1(1, aNPnt2-1);
289 math_Vector f2(1, aNPnt2-1);
291 Standard_Real aXm = 0.5*(theUpper + theLower);
292 Standard_Real aXr = 0.5*(theUpper - theLower);
294 // Compute Gauss quadrature
298 for (i = 2; i < aNPnt2; i += 2) {
299 aDx = aXr*theKronrodP.Value(i);
301 if (!theFunction.Value(aXm + aDx, aVal1) ||
302 !theFunction.Value(aXm - aDx, aVal2)) {
303 IsDone = Standard_False;
309 aGaussVal += (aVal1 + aVal2)*theGaussW.Value(i/2);
310 theValue += (aVal1 + aVal2)*theKronrodW.Value(i);
313 // Compute value in the middle point.
314 if (!theFunction.Value(aXm, aVal1)) {
315 IsDone = Standard_False;
319 Standard_Real fc = aVal1;
320 theValue += aVal1*theKronrodW.Value(aNPnt2);
323 aGaussVal += aVal1*theGaussW.Value(aNPnt2/2);
325 // Compute Kronrod quadrature
326 for (i = 1; i < aNPnt2; i += 2) {
327 aDx = aXr*theKronrodP.Value(i);
329 if (!theFunction.Value(aXm + aDx, aVal1) ||
330 !theFunction.Value(aXm - aDx, aVal2)) {
331 IsDone = Standard_False;
337 theValue += (aVal1 + aVal2)*theKronrodW.Value(i);
340 Standard_Real mean = 0.5*theValue;
342 Standard_Real asc = Abs(fc-mean)*theKronrodW.Value(aNPnt2);
343 for(i = 1; i < aNPnt2; ++i) {
344 asc += theKronrodW.Value(i)*(Abs(f1(i) - mean) + Abs(f2(i) - mean));
351 // Compute the error and the new number of Kronrod points.
353 theError = Abs(theValue - aGaussVal);
355 Standard_Real scale =1.;
356 if(asc != 0. && theError != 0.) scale = Pow((200.*theError/asc), 1.5);
357 if(scale < 1.) theError = Min(theError, asc*scale);
359 //theFunction.GetStateNumber();
361 return Standard_True;