src/math/math_KronrodSingleIntegration.cxx

   1 // Created on: 2005-12-08
   2 // Created by: Sergey KHROMOV
   3 // Copyright (c) 2005-2012 OPEN CASCADE SAS
   4 //
   5 // The content of this file is subject to the Open CASCADE Technology Public
   6 // License Version 6.5 (the "License"). You may not use the content of this file
   7 // except in compliance with the License. Please obtain a copy of the License
   8 // at http://www.opencascade.org and read it completely before using this file.
   9 //
  10 // The Initial Developer of the Original Code is Open CASCADE S.A.S., having its
  11 // main offices at: 1, place des Freres Montgolfier, 78280 Guyancourt, France.
  12 //
  13 // The Original Code and all software distributed under the License is
  14 // distributed on an "AS IS" basis, without warranty of any kind, and the
  15 // Initial Developer hereby disclaims all such warranties, including without
  16 // limitation, any warranties of merchantability, fitness for a particular
  17 // purpose or non-infringement. Please see the License for the specific terms
  18 // and conditions governing the rights and limitations under the License.
  19
  20
  21
  22 #include <math_KronrodSingleIntegration.ixx>
  23 #include <math_Vector.hxx>
  24 #include <math.hxx>
  25 #include <TColStd_SequenceOfReal.hxx>
  26
  27
  28 //==========================================================================
  29 //function : An empty constructor.
  30 //==========================================================================
  31
  32 math_KronrodSingleIntegration::math_KronrodSingleIntegration() :
  33        myIsDone(Standard_False),
  34        myValue(0.),
  35        myErrorReached(0.),
  36        myNbPntsReached(0),
  37        myNbIterReached(0)
  38 {
  39 }
  40
  41 //==========================================================================
  42 //function : Constructor
  43 //
  44 //==========================================================================
  45
  46 math_KronrodSingleIntegration::math_KronrodSingleIntegration
  47                               (      math_Function    &theFunction,
  48                                const Standard_Real     theLower,
  49                                const Standard_Real     theUpper,
  50                                const Standard_Integer  theNbPnts):
  51   myIsDone(Standard_False),
  52   myValue(0.),
  53   myErrorReached(0.),
  54   myNbPntsReached(0)
  55 {
  56   Perform(theFunction, theLower, theUpper, theNbPnts);
  57 }
  58
  59 //==========================================================================
  60 //function : Constructor
  61 //
  62 //==========================================================================
  63
  64 math_KronrodSingleIntegration::math_KronrodSingleIntegration
  65                               (      math_Function    &theFunction,
  66                                const Standard_Real     theLower,
  67                                const Standard_Real     theUpper,
  68                                const Standard_Integer  theNbPnts,
  69                                const Standard_Real     theTolerance,
  70                                const Standard_Integer  theMaxNbIter) :
  71   myIsDone(Standard_False),
  72   myValue(0.),
  73   myErrorReached(0.),
  74   myNbPntsReached(0)
  75 {
  76   Perform(theFunction, theLower, theUpper, theNbPnts,
  77           theTolerance, theMaxNbIter);
  78 }
  79
  80 //==========================================================================
  81 //function : Perform
  82 //           Computation of the integral.
  83 //==========================================================================
  84
  85 void math_KronrodSingleIntegration::Perform
  86                               (      math_Function    &theFunction,
  87                                const Standard_Real     theLower,
  88                                const Standard_Real     theUpper,
  89                                const Standard_Integer  theNbPnts)
  90 {
  91   //const Standard_Real aMinVol = Epsilon(1.);
  92   const Standard_Real aPtol = 1.e-9;
  93   myNbIterReached = 0;
  94
  95   if (theNbPnts    <  3 ) {
  96     myIsDone = Standard_False;
  97     return;
  98   }
  99
 100   if(theUpper - theLower < aPtol) {
 101     myIsDone = Standard_False;
 102     return;
 103   }
 104
 105   // Get an odd value of number of initial points.
 106   myNbPntsReached = (theNbPnts%2 == 0) ? theNbPnts + 1 : theNbPnts;
 107   myErrorReached  = RealLast();
 108
 109   Standard_Integer aNGauss = myNbPntsReached/2;
 110   math_Vector      aKronrodP(1, myNbPntsReached);
 111   math_Vector      aKronrodW(1, myNbPntsReached);
 112   math_Vector      aGaussP(1, aNGauss);
 113   math_Vector      aGaussW(1, aNGauss);
 114
 115
 116   if (!math::KronrodPointsAndWeights(myNbPntsReached, aKronrodP, aKronrodW) ||
 117       !math::OrderedGaussPointsAndWeights(aNGauss, aGaussP, aGaussW)) {
 118     myIsDone = Standard_False;
 119     return;
 120   }
 121
 122   myIsDone = GKRule(theFunction, theLower, theUpper, aGaussP, aGaussW, aKronrodP, aKronrodW,
 123                     myValue, myErrorReached);
 124
 125   if(!myIsDone) return;
 126
 127   //Standard_Real anAbsVal = Abs(myValue);
 128
 129   myAbsolutError = myErrorReached;
 130
 131   //if (anAbsVal > aMinVol)
 132     //myErrorReached /= anAbsVal;
 133
 134   myNbIterReached++;
 135
 136 }
 137
 138 //=======================================================================
 139 //function : Perform
 140
 141 //purpose  :
 142 //=======================================================================
 143
 144 void math_KronrodSingleIntegration::Perform
 145                               (      math_Function    &theFunction,
 146                                const Standard_Real     theLower,
 147                                const Standard_Real     theUpper,
 148                                const Standard_Integer  theNbPnts,
 149                                const Standard_Real     theTolerance,
 150                                const Standard_Integer  theMaxNbIter)
 151 {
 152   Standard_Real aMinVol = Epsilon(1.);
 153   myNbIterReached = 0;
 154
 155   // Check prerequisites.
 156   if (theNbPnts    <  3 ||
 157       theTolerance <= 0.) {
 158     myIsDone = Standard_False;
 159     return;
 160   }
 161   // Get an odd value of number of initial points.
 162   myNbPntsReached = (theNbPnts%2 == 0) ? theNbPnts + 1 : theNbPnts;
 163
 164   Standard_Integer aNGauss = myNbPntsReached/2;
 165   math_Vector      aKronrodP(1, myNbPntsReached);
 166   math_Vector      aKronrodW(1, myNbPntsReached);
 167   math_Vector      aGaussP(1, aNGauss);
 168   math_Vector      aGaussW(1, aNGauss);
 169
 170   if (!math::KronrodPointsAndWeights(myNbPntsReached, aKronrodP, aKronrodW) ||
 171       !math::OrderedGaussPointsAndWeights(aNGauss, aGaussP, aGaussW)) {
 172     myIsDone = Standard_False;
 173     return;
 174   }
 175
 176   //First iteration
 177   myIsDone = GKRule(theFunction, theLower, theUpper, aGaussP, aGaussW, aKronrodP, aKronrodW,
 178                     myValue, myErrorReached);
 179
 180   if(!myIsDone) return;
 181
 182   Standard_Real anAbsVal = Abs(myValue);
 183
 184   myAbsolutError = myErrorReached;
 185   if (anAbsVal > aMinVol)
 186     myErrorReached /= anAbsVal;
 187
 188   myNbIterReached++;
 189
 190   if(myErrorReached <= theTolerance) return;
 191   if(myNbIterReached >= theMaxNbIter) return;
 192
 193   TColStd_SequenceOfReal anIntervals;
 194   TColStd_SequenceOfReal anErrors;
 195   TColStd_SequenceOfReal aValues;
 196
 197   anIntervals.Append(theLower);
 198   anIntervals.Append(theUpper);
 199
 200
 201   anErrors.Append(myAbsolutError);
 202   aValues.Append(myValue);
 203
 204   Standard_Integer i, nint, nbints;
 205
 206   Standard_Real maxerr;
 207   Standard_Integer count = 0;
 208
 209   while(myErrorReached > theTolerance && myNbIterReached < theMaxNbIter) {
 210     //Searching interval with max error
 211     nbints = anIntervals.Length() - 1;
 212     nint = 0;
 213     maxerr = 0.;
 214     for(i = 1; i <= nbints; ++i) {
 215       if(anErrors(i) > maxerr) {
 216         maxerr = anErrors(i);
 217         nint = i;
 218       }
 219     }
 220
 221     Standard_Real a = anIntervals(nint);
 222     Standard_Real b = anIntervals(nint+1);
 223     Standard_Real c = 0.5*(a + b);
 224
 225     Standard_Real v1, v2, e1, e2;
 226
 227     myIsDone = GKRule(theFunction, a, c, aGaussP, aGaussW, aKronrodP, aKronrodW,
 228                     v1, e1);
 229
 230     if(!myIsDone) return;
 231
 232     myIsDone = GKRule(theFunction, c, b, aGaussP, aGaussW, aKronrodP, aKronrodW,
 233                     v2, e2);
 234
 235     if(!myIsDone) return;
 236
 237     myNbIterReached++;
 238
 239     Standard_Real deltav = v1 + v2 - aValues(nint);
 240     myValue += deltav;
 241     if(Abs(deltav) <= Epsilon(Abs(myValue))) ++count;
 242
 243     Standard_Real deltae = e1 + e2 - anErrors(nint);
 244     myAbsolutError +=  deltae;
 245     if(myAbsolutError <= Epsilon(Abs(myValue))) ++count;
 246
 247     if(Abs(myValue) > aMinVol) myErrorReached = myAbsolutError/Abs(myValue);
 248     else myErrorReached = myAbsolutError;
 249
 250
 251     if(count > 50) return;
 252
 253     //Inserting new interval
 254
 255     anIntervals.InsertAfter(nint, c);
 256     anErrors(nint) = e1;
 257     anErrors.InsertAfter(nint, e2);
 258     aValues(nint) = v1;
 259     aValues.InsertAfter(nint, v2);
 260
 261   }
 262
 263 }
 264
 265 //=======================================================================
 266 //function : GKRule
 267 //purpose  :
 268 //=======================================================================
 269
 270 Standard_Boolean math_KronrodSingleIntegration::GKRule(
 271                                      math_Function    &theFunction,
 272                                const Standard_Real     theLower,
 273                                const Standard_Real     theUpper,
 274                                const math_Vector&      /*theGaussP*/,
 275                                const math_Vector&      theGaussW,
 276                                const math_Vector&      theKronrodP,
 277                                const math_Vector&      theKronrodW,
 278                                      Standard_Real&    theValue,
 279                                      Standard_Real&    theError)
 280 {
 281
 282   Standard_Boolean IsDone;
 283
 284   Standard_Integer aNKronrod = theKronrodP.Length();
 285
 286   Standard_Real    aGaussVal;
 287   Standard_Integer aNPnt2 = (aNKronrod + 1)/2;
 288   Standard_Integer i;
 289   Standard_Real    aDx;
 290   Standard_Real    aVal1;
 291   Standard_Real    aVal2;
 292
 293   math_Vector      f1(1, aNPnt2-1);
 294   math_Vector      f2(1, aNPnt2-1);
 295
 296   Standard_Real    aXm = 0.5*(theUpper + theLower);
 297   Standard_Real    aXr = 0.5*(theUpper - theLower);
 298
 299   // Compute Gauss quadrature
 300   aGaussVal = 0.;
 301   theValue   = 0.;
 302
 303   for (i = 2; i < aNPnt2; i += 2) {
 304     aDx = aXr*theKronrodP.Value(i);
 305
 306     if (!theFunction.Value(aXm + aDx, aVal1) ||
 307         !theFunction.Value(aXm - aDx, aVal2)) {
 308       IsDone = Standard_False;
 309       return IsDone;
 310     }
 311
 312     f1(i) = aVal1;
 313     f2(i) = aVal2;
 314     aGaussVal += (aVal1 + aVal2)*theGaussW.Value(i/2);
 315     theValue   += (aVal1 + aVal2)*theKronrodW.Value(i);
 316   }
 317
 318   // Compute value in the middle point.
 319   if (!theFunction.Value(aXm, aVal1)) {
 320     IsDone = Standard_False;
 321     return IsDone;
 322   }
 323
 324   Standard_Real fc = aVal1;
 325   theValue += aVal1*theKronrodW.Value(aNPnt2);
 326
 327   if (i == aNPnt2)
 328     aGaussVal += aVal1*theGaussW.Value(aNPnt2/2);
 329
 330   // Compute Kronrod quadrature
 331   for (i = 1; i < aNPnt2; i += 2) {
 332     aDx = aXr*theKronrodP.Value(i);
 333
 334     if (!theFunction.Value(aXm + aDx, aVal1) ||
 335         !theFunction.Value(aXm - aDx, aVal2)) {
 336       IsDone = Standard_False;
 337       return IsDone;
 338     }
 339
 340     f1(i) = aVal1;
 341     f2(i) = aVal2;
 342     theValue += (aVal1 + aVal2)*theKronrodW.Value(i);
 343   }
 344
 345   Standard_Real mean = 0.5*theValue;
 346
 347   Standard_Real asc = Abs(fc-mean)*theKronrodW.Value(aNPnt2);
 348   for(i = 1; i < aNPnt2; ++i) {
 349     asc += theKronrodW.Value(i)*(Abs(f1(i) - mean) + Abs(f2(i) - mean));
 350   }
 351
 352   asc *= aXr;
 353   theValue   *= aXr;
 354   aGaussVal *= aXr;
 355
 356   // Compute the error and the new number of Kronrod points.
 357
 358   theError = Abs(theValue - aGaussVal);
 359
 360   Standard_Real scale =1.;
 361   if(asc != 0. && theError != 0.) scale = Pow((200.*theError/asc), 1.5);
 362   if(scale < 1.) theError = Min(theError, asc*scale);
 363
 364   //theFunction.GetStateNumber();
 365
 366   return Standard_True;
 367
 368 }
 369