src/math/math_KronrodSingleIntegration.cxx

   1 // Created on: 2005-12-08
   2 // Created by: Sergey KHROMOV
   3 // Copyright (c) 2005-2014 OPEN CASCADE SAS
   4 //
   5 // This file is part of Open CASCADE Technology software library.
   6 //
   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.
  12 //
  13 // Alternatively, this file may be used under the terms of Open CASCADE
  14 // commercial license or contractual agreement.
  15
  16
  17 #include <math.hxx>
  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>
  23
  24 //==========================================================================
  25 //function : An empty constructor.
  26 //==========================================================================
  27 math_KronrodSingleIntegration::math_KronrodSingleIntegration() :
  28        myIsDone(Standard_False),
  29        myValue(0.),
  30        myErrorReached(0.),
  31        myNbPntsReached(0),
  32        myNbIterReached(0)
  33 {
  34 }
  35
  36 //==========================================================================
  37 //function : Constructor
  38 //
  39 //==========================================================================
  40
  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),
  47   myValue(0.),
  48   myErrorReached(0.),
  49   myNbPntsReached(0)
  50 {
  51   Perform(theFunction, theLower, theUpper, theNbPnts);
  52 }
  53
  54 //==========================================================================
  55 //function : Constructor
  56 //
  57 //==========================================================================
  58
  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),
  67   myValue(0.),
  68   myErrorReached(0.),
  69   myNbPntsReached(0)
  70 {
  71   Perform(theFunction, theLower, theUpper, theNbPnts,
  72           theTolerance, theMaxNbIter);
  73 }
  74
  75 //==========================================================================
  76 //function : Perform
  77 //           Computation of the integral.
  78 //==========================================================================
  79
  80 void math_KronrodSingleIntegration::Perform
  81                               (      math_Function    &theFunction,
  82                                const Standard_Real     theLower,
  83                                const Standard_Real     theUpper,
  84                                const Standard_Integer  theNbPnts)
  85 {
  86   //const Standard_Real aMinVol = Epsilon(1.);
  87   const Standard_Real aPtol = 1.e-9;
  88   myNbIterReached = 0;
  89
  90   if (theNbPnts    <  3 ) {
  91     myIsDone = Standard_False;
  92     return;
  93   }
  94
  95   if(theUpper - theLower < aPtol) {
  96     myIsDone = Standard_False;
  97     return;
  98   }
  99
 100   // Get an odd value of number of initial points.
 101   myNbPntsReached = (theNbPnts%2 == 0) ? theNbPnts + 1 : theNbPnts;
 102   myErrorReached  = RealLast();
 103
 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);
 109
 110
 111   if (!math::KronrodPointsAndWeights(myNbPntsReached, aKronrodP, aKronrodW) ||
 112       !math::OrderedGaussPointsAndWeights(aNGauss, aGaussP, aGaussW)) {
 113     myIsDone = Standard_False;
 114     return;
 115   }
 116
 117   myIsDone = GKRule(theFunction, theLower, theUpper, aGaussP, aGaussW, aKronrodP, aKronrodW,
 118                     myValue, myErrorReached);
 119
 120   if(!myIsDone) return;
 121
 122   //Standard_Real anAbsVal = Abs(myValue);
 123
 124   myAbsolutError = myErrorReached;
 125
 126   //if (anAbsVal > aMinVol)
 127     //myErrorReached /= anAbsVal;
 128
 129   myNbIterReached++;
 130
 131 }
 132
 133 //=======================================================================
 134 //function : Perform
 135
 136 //purpose  :
 137 //=======================================================================
 138
 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)
 146 {
 147   Standard_Real aMinVol = Epsilon(1.);
 148   myNbIterReached = 0;
 149
 150   // Check prerequisites.
 151   if (theNbPnts    <  3 ||
 152       theTolerance <= 0.) {
 153     myIsDone = Standard_False;
 154     return;
 155   }
 156   // Get an odd value of number of initial points.
 157   myNbPntsReached = (theNbPnts%2 == 0) ? theNbPnts + 1 : theNbPnts;
 158
 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);
 164
 165   if (!math::KronrodPointsAndWeights(myNbPntsReached, aKronrodP, aKronrodW) ||
 166       !math::OrderedGaussPointsAndWeights(aNGauss, aGaussP, aGaussW)) {
 167     myIsDone = Standard_False;
 168     return;
 169   }
 170
 171   //First iteration
 172   myIsDone = GKRule(theFunction, theLower, theUpper, aGaussP, aGaussW, aKronrodP, aKronrodW,
 173                     myValue, myErrorReached);
 174
 175   if(!myIsDone) return;
 176
 177   Standard_Real anAbsVal = Abs(myValue);
 178
 179   myAbsolutError = myErrorReached;
 180   if (anAbsVal > aMinVol)
 181     myErrorReached /= anAbsVal;
 182
 183   myNbIterReached++;
 184
 185   if(myErrorReached <= theTolerance) return;
 186   if(myNbIterReached >= theMaxNbIter) return;
 187
 188   TColStd_SequenceOfReal anIntervals;
 189   TColStd_SequenceOfReal anErrors;
 190   TColStd_SequenceOfReal aValues;
 191
 192   anIntervals.Append(theLower);
 193   anIntervals.Append(theUpper);
 194
 195
 196   anErrors.Append(myAbsolutError);
 197   aValues.Append(myValue);
 198
 199   Standard_Integer i, nint, nbints;
 200
 201   Standard_Real maxerr;
 202   Standard_Integer count = 0;
 203
 204   while(myErrorReached > theTolerance && myNbIterReached < theMaxNbIter) {
 205     //Searching interval with max error
 206     nbints = anIntervals.Length() - 1;
 207     nint = 0;
 208     maxerr = 0.;
 209     for(i = 1; i <= nbints; ++i) {
 210       if(anErrors(i) > maxerr) {
 211         maxerr = anErrors(i);
 212         nint = i;
 213       }
 214     }
 215
 216     Standard_Real a = anIntervals(nint);
 217     Standard_Real b = anIntervals(nint+1);
 218     Standard_Real c = 0.5*(a + b);
 219
 220     Standard_Real v1, v2, e1, e2;
 221
 222     myIsDone = GKRule(theFunction, a, c, aGaussP, aGaussW, aKronrodP, aKronrodW,
 223                     v1, e1);
 224
 225     if(!myIsDone) return;
 226
 227     myIsDone = GKRule(theFunction, c, b, aGaussP, aGaussW, aKronrodP, aKronrodW,
 228                     v2, e2);
 229
 230     if(!myIsDone) return;
 231
 232     myNbIterReached++;
 233
 234     Standard_Real deltav = v1 + v2 - aValues(nint);
 235     myValue += deltav;
 236     if(Abs(deltav) <= Epsilon(Abs(myValue))) ++count;
 237
 238     Standard_Real deltae = e1 + e2 - anErrors(nint);
 239     myAbsolutError +=  deltae;
 240     if(myAbsolutError <= Epsilon(Abs(myValue))) ++count;
 241
 242     if(Abs(myValue) > aMinVol) myErrorReached = myAbsolutError/Abs(myValue);
 243     else myErrorReached = myAbsolutError;
 244
 245
 246     if(count > 50) return;
 247
 248     //Inserting new interval
 249
 250     anIntervals.InsertAfter(nint, c);
 251     anErrors(nint) = e1;
 252     anErrors.InsertAfter(nint, e2);
 253     aValues(nint) = v1;
 254     aValues.InsertAfter(nint, v2);
 255
 256   }
 257
 258 }
 259
 260 //=======================================================================
 261 //function : GKRule
 262 //purpose  :
 263 //=======================================================================
 264
 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)
 275 {
 276
 277   Standard_Boolean IsDone;
 278
 279   Standard_Integer aNKronrod = theKronrodP.Length();
 280
 281   Standard_Real    aGaussVal;
 282   Standard_Integer aNPnt2 = (aNKronrod + 1)/2;
 283   Standard_Integer i;
 284   Standard_Real    aDx;
 285   Standard_Real    aVal1;
 286   Standard_Real    aVal2;
 287
 288   math_Vector      f1(1, aNPnt2-1);
 289   math_Vector      f2(1, aNPnt2-1);
 290
 291   Standard_Real    aXm = 0.5*(theUpper + theLower);
 292   Standard_Real    aXr = 0.5*(theUpper - theLower);
 293
 294   // Compute Gauss quadrature
 295   aGaussVal = 0.;
 296   theValue   = 0.;
 297
 298   for (i = 2; i < aNPnt2; i += 2) {
 299     aDx = aXr*theKronrodP.Value(i);
 300
 301     if (!theFunction.Value(aXm + aDx, aVal1) ||
 302         !theFunction.Value(aXm - aDx, aVal2)) {
 303       IsDone = Standard_False;
 304       return IsDone;
 305     }
 306
 307     f1(i) = aVal1;
 308     f2(i) = aVal2;
 309     aGaussVal += (aVal1 + aVal2)*theGaussW.Value(i/2);
 310     theValue   += (aVal1 + aVal2)*theKronrodW.Value(i);
 311   }
 312
 313   // Compute value in the middle point.
 314   if (!theFunction.Value(aXm, aVal1)) {
 315     IsDone = Standard_False;
 316     return IsDone;
 317   }
 318
 319   Standard_Real fc = aVal1;
 320   theValue += aVal1*theKronrodW.Value(aNPnt2);
 321
 322   if (i == aNPnt2)
 323     aGaussVal += aVal1*theGaussW.Value(aNPnt2/2);
 324
 325   // Compute Kronrod quadrature
 326   for (i = 1; i < aNPnt2; i += 2) {
 327     aDx = aXr*theKronrodP.Value(i);
 328
 329     if (!theFunction.Value(aXm + aDx, aVal1) ||
 330         !theFunction.Value(aXm - aDx, aVal2)) {
 331       IsDone = Standard_False;
 332       return IsDone;
 333     }
 334
 335     f1(i) = aVal1;
 336     f2(i) = aVal2;
 337     theValue += (aVal1 + aVal2)*theKronrodW.Value(i);
 338   }
 339
 340   Standard_Real mean = 0.5*theValue;
 341
 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));
 345   }
 346
 347   asc *= aXr;
 348   theValue   *= aXr;
 349   aGaussVal *= aXr;
 350
 351   // Compute the error and the new number of Kronrod points.
 352
 353   theError = Abs(theValue - aGaussVal);
 354
 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);
 358
 359   //theFunction.GetStateNumber();
 360
 361   return Standard_True;
 362
 363 }
 364