mirror of
https://git.dev.opencascade.org/repos/occt.git
synced 2025-04-08 18:40:55 +03:00
b9280b8b27
Removed unused exception classes OSD_Exception_FLT_DIVIDE_BY_ZERO, OSD_Exception_INT_DIVIDE_BY_ZERO, OSD_Exception_FLT_DENORMAL_OPERAND, OSD_Exception_FLT_INEXACT_RESULT, OSD_Exception_FLT_INVALID_OPERATION, OSD_Exception_FLT_OVERFLOW, OSD_Exception_FLT_STACK_CHECK, OSD_Exception_FLT_UNDERFLOW.
364 lines
10 KiB
C++
364 lines
10 KiB
C++
// Created on: 2005-12-08
|
|
// Created by: Sergey KHROMOV
|
|
// Copyright (c) 2005-2014 OPEN CASCADE SAS
|
|
//
|
|
// This file is part of Open CASCADE Technology software library.
|
|
//
|
|
// This library is free software; you can redistribute it and/or modify it under
|
|
// the terms of the GNU Lesser General Public License version 2.1 as published
|
|
// by the Free Software Foundation, with special exception defined in the file
|
|
// OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT
|
|
// distribution for complete text of the license and disclaimer of any warranty.
|
|
//
|
|
// Alternatively, this file may be used under the terms of Open CASCADE
|
|
// commercial license or contractual agreement.
|
|
|
|
|
|
#include <math.hxx>
|
|
#include <math_Function.hxx>
|
|
#include <math_KronrodSingleIntegration.hxx>
|
|
#include <math_Vector.hxx>
|
|
#include <TColStd_SequenceOfReal.hxx>
|
|
|
|
//==========================================================================
|
|
//function : An empty constructor.
|
|
//==========================================================================
|
|
math_KronrodSingleIntegration::math_KronrodSingleIntegration() :
|
|
myIsDone(Standard_False),
|
|
myValue(0.),
|
|
myErrorReached(0.),
|
|
myNbPntsReached(0),
|
|
myNbIterReached(0)
|
|
{
|
|
}
|
|
|
|
//==========================================================================
|
|
//function : Constructor
|
|
//
|
|
//==========================================================================
|
|
|
|
math_KronrodSingleIntegration::math_KronrodSingleIntegration
|
|
( math_Function &theFunction,
|
|
const Standard_Real theLower,
|
|
const Standard_Real theUpper,
|
|
const Standard_Integer theNbPnts):
|
|
myIsDone(Standard_False),
|
|
myValue(0.),
|
|
myErrorReached(0.),
|
|
myNbPntsReached(0)
|
|
{
|
|
Perform(theFunction, theLower, theUpper, theNbPnts);
|
|
}
|
|
|
|
//==========================================================================
|
|
//function : Constructor
|
|
//
|
|
//==========================================================================
|
|
|
|
math_KronrodSingleIntegration::math_KronrodSingleIntegration
|
|
( math_Function &theFunction,
|
|
const Standard_Real theLower,
|
|
const Standard_Real theUpper,
|
|
const Standard_Integer theNbPnts,
|
|
const Standard_Real theTolerance,
|
|
const Standard_Integer theMaxNbIter) :
|
|
myIsDone(Standard_False),
|
|
myValue(0.),
|
|
myErrorReached(0.),
|
|
myNbPntsReached(0)
|
|
{
|
|
Perform(theFunction, theLower, theUpper, theNbPnts,
|
|
theTolerance, theMaxNbIter);
|
|
}
|
|
|
|
//==========================================================================
|
|
//function : Perform
|
|
// Computation of the integral.
|
|
//==========================================================================
|
|
|
|
void math_KronrodSingleIntegration::Perform
|
|
( math_Function &theFunction,
|
|
const Standard_Real theLower,
|
|
const Standard_Real theUpper,
|
|
const Standard_Integer theNbPnts)
|
|
{
|
|
//const Standard_Real aMinVol = Epsilon(1.);
|
|
const Standard_Real aPtol = 1.e-9;
|
|
myNbIterReached = 0;
|
|
|
|
if (theNbPnts < 3 ) {
|
|
myIsDone = Standard_False;
|
|
return;
|
|
}
|
|
|
|
if(theUpper - theLower < aPtol) {
|
|
myIsDone = Standard_False;
|
|
return;
|
|
}
|
|
|
|
// Get an odd value of number of initial points.
|
|
myNbPntsReached = (theNbPnts%2 == 0) ? theNbPnts + 1 : theNbPnts;
|
|
myErrorReached = RealLast();
|
|
|
|
Standard_Integer aNGauss = myNbPntsReached/2;
|
|
math_Vector aKronrodP(1, myNbPntsReached);
|
|
math_Vector aKronrodW(1, myNbPntsReached);
|
|
math_Vector aGaussP(1, aNGauss);
|
|
math_Vector aGaussW(1, aNGauss);
|
|
|
|
|
|
if (!math::KronrodPointsAndWeights(myNbPntsReached, aKronrodP, aKronrodW) ||
|
|
!math::OrderedGaussPointsAndWeights(aNGauss, aGaussP, aGaussW)) {
|
|
myIsDone = Standard_False;
|
|
return;
|
|
}
|
|
|
|
myIsDone = GKRule(theFunction, theLower, theUpper, aGaussP, aGaussW, aKronrodP, aKronrodW,
|
|
myValue, myErrorReached);
|
|
|
|
if(!myIsDone) return;
|
|
|
|
//Standard_Real anAbsVal = Abs(myValue);
|
|
|
|
myAbsolutError = myErrorReached;
|
|
|
|
//if (anAbsVal > aMinVol)
|
|
//myErrorReached /= anAbsVal;
|
|
|
|
myNbIterReached++;
|
|
|
|
}
|
|
|
|
//=======================================================================
|
|
//function : Perform
|
|
|
|
//purpose :
|
|
//=======================================================================
|
|
|
|
void math_KronrodSingleIntegration::Perform
|
|
( math_Function &theFunction,
|
|
const Standard_Real theLower,
|
|
const Standard_Real theUpper,
|
|
const Standard_Integer theNbPnts,
|
|
const Standard_Real theTolerance,
|
|
const Standard_Integer theMaxNbIter)
|
|
{
|
|
Standard_Real aMinVol = Epsilon(1.);
|
|
myNbIterReached = 0;
|
|
|
|
// Check prerequisites.
|
|
if (theNbPnts < 3 ||
|
|
theTolerance <= 0.) {
|
|
myIsDone = Standard_False;
|
|
return;
|
|
}
|
|
// Get an odd value of number of initial points.
|
|
myNbPntsReached = (theNbPnts%2 == 0) ? theNbPnts + 1 : theNbPnts;
|
|
|
|
Standard_Integer aNGauss = myNbPntsReached/2;
|
|
math_Vector aKronrodP(1, myNbPntsReached);
|
|
math_Vector aKronrodW(1, myNbPntsReached);
|
|
math_Vector aGaussP(1, aNGauss);
|
|
math_Vector aGaussW(1, aNGauss);
|
|
|
|
if (!math::KronrodPointsAndWeights(myNbPntsReached, aKronrodP, aKronrodW) ||
|
|
!math::OrderedGaussPointsAndWeights(aNGauss, aGaussP, aGaussW)) {
|
|
myIsDone = Standard_False;
|
|
return;
|
|
}
|
|
|
|
//First iteration
|
|
myIsDone = GKRule(theFunction, theLower, theUpper, aGaussP, aGaussW, aKronrodP, aKronrodW,
|
|
myValue, myErrorReached);
|
|
|
|
if(!myIsDone) return;
|
|
|
|
Standard_Real anAbsVal = Abs(myValue);
|
|
|
|
myAbsolutError = myErrorReached;
|
|
if (anAbsVal > aMinVol)
|
|
myErrorReached /= anAbsVal;
|
|
|
|
myNbIterReached++;
|
|
|
|
if(myErrorReached <= theTolerance) return;
|
|
if(myNbIterReached >= theMaxNbIter) return;
|
|
|
|
TColStd_SequenceOfReal anIntervals;
|
|
TColStd_SequenceOfReal anErrors;
|
|
TColStd_SequenceOfReal aValues;
|
|
|
|
anIntervals.Append(theLower);
|
|
anIntervals.Append(theUpper);
|
|
|
|
|
|
anErrors.Append(myAbsolutError);
|
|
aValues.Append(myValue);
|
|
|
|
Standard_Integer i, nint, nbints;
|
|
|
|
Standard_Real maxerr;
|
|
Standard_Integer count = 0;
|
|
|
|
while(myErrorReached > theTolerance && myNbIterReached < theMaxNbIter) {
|
|
//Searching interval with max error
|
|
nbints = anIntervals.Length() - 1;
|
|
nint = 0;
|
|
maxerr = 0.;
|
|
for(i = 1; i <= nbints; ++i) {
|
|
if(anErrors(i) > maxerr) {
|
|
maxerr = anErrors(i);
|
|
nint = i;
|
|
}
|
|
}
|
|
|
|
Standard_Real a = anIntervals(nint);
|
|
Standard_Real b = anIntervals(nint+1);
|
|
Standard_Real c = 0.5*(a + b);
|
|
|
|
Standard_Real v1, v2, e1, e2;
|
|
|
|
myIsDone = GKRule(theFunction, a, c, aGaussP, aGaussW, aKronrodP, aKronrodW,
|
|
v1, e1);
|
|
|
|
if(!myIsDone) return;
|
|
|
|
myIsDone = GKRule(theFunction, c, b, aGaussP, aGaussW, aKronrodP, aKronrodW,
|
|
v2, e2);
|
|
|
|
if(!myIsDone) return;
|
|
|
|
myNbIterReached++;
|
|
|
|
Standard_Real deltav = v1 + v2 - aValues(nint);
|
|
myValue += deltav;
|
|
if(Abs(deltav) <= Epsilon(Abs(myValue))) ++count;
|
|
|
|
Standard_Real deltae = e1 + e2 - anErrors(nint);
|
|
myAbsolutError += deltae;
|
|
if(myAbsolutError <= Epsilon(Abs(myValue))) ++count;
|
|
|
|
if(Abs(myValue) > aMinVol) myErrorReached = myAbsolutError/Abs(myValue);
|
|
else myErrorReached = myAbsolutError;
|
|
|
|
|
|
if(count > 50) return;
|
|
|
|
//Inserting new interval
|
|
|
|
anIntervals.InsertAfter(nint, c);
|
|
anErrors(nint) = e1;
|
|
anErrors.InsertAfter(nint, e2);
|
|
aValues(nint) = v1;
|
|
aValues.InsertAfter(nint, v2);
|
|
|
|
}
|
|
|
|
}
|
|
|
|
//=======================================================================
|
|
//function : GKRule
|
|
//purpose :
|
|
//=======================================================================
|
|
|
|
Standard_Boolean math_KronrodSingleIntegration::GKRule(
|
|
math_Function &theFunction,
|
|
const Standard_Real theLower,
|
|
const Standard_Real theUpper,
|
|
const math_Vector& /*theGaussP*/,
|
|
const math_Vector& theGaussW,
|
|
const math_Vector& theKronrodP,
|
|
const math_Vector& theKronrodW,
|
|
Standard_Real& theValue,
|
|
Standard_Real& theError)
|
|
{
|
|
|
|
Standard_Boolean IsDone;
|
|
|
|
Standard_Integer aNKronrod = theKronrodP.Length();
|
|
|
|
Standard_Real aGaussVal;
|
|
Standard_Integer aNPnt2 = (aNKronrod + 1)/2;
|
|
Standard_Integer i;
|
|
Standard_Real aDx;
|
|
Standard_Real aVal1;
|
|
Standard_Real aVal2;
|
|
|
|
math_Vector f1(1, aNPnt2-1);
|
|
math_Vector f2(1, aNPnt2-1);
|
|
|
|
Standard_Real aXm = 0.5*(theUpper + theLower);
|
|
Standard_Real aXr = 0.5*(theUpper - theLower);
|
|
|
|
// Compute Gauss quadrature
|
|
aGaussVal = 0.;
|
|
theValue = 0.;
|
|
|
|
for (i = 2; i < aNPnt2; i += 2) {
|
|
aDx = aXr*theKronrodP.Value(i);
|
|
|
|
if (!theFunction.Value(aXm + aDx, aVal1) ||
|
|
!theFunction.Value(aXm - aDx, aVal2)) {
|
|
IsDone = Standard_False;
|
|
return IsDone;
|
|
}
|
|
|
|
f1(i) = aVal1;
|
|
f2(i) = aVal2;
|
|
aGaussVal += (aVal1 + aVal2)*theGaussW.Value(i/2);
|
|
theValue += (aVal1 + aVal2)*theKronrodW.Value(i);
|
|
}
|
|
|
|
// Compute value in the middle point.
|
|
if (!theFunction.Value(aXm, aVal1)) {
|
|
IsDone = Standard_False;
|
|
return IsDone;
|
|
}
|
|
|
|
Standard_Real fc = aVal1;
|
|
theValue += aVal1*theKronrodW.Value(aNPnt2);
|
|
|
|
if (i == aNPnt2)
|
|
aGaussVal += aVal1*theGaussW.Value(aNPnt2/2);
|
|
|
|
// Compute Kronrod quadrature
|
|
for (i = 1; i < aNPnt2; i += 2) {
|
|
aDx = aXr*theKronrodP.Value(i);
|
|
|
|
if (!theFunction.Value(aXm + aDx, aVal1) ||
|
|
!theFunction.Value(aXm - aDx, aVal2)) {
|
|
IsDone = Standard_False;
|
|
return IsDone;
|
|
}
|
|
|
|
f1(i) = aVal1;
|
|
f2(i) = aVal2;
|
|
theValue += (aVal1 + aVal2)*theKronrodW.Value(i);
|
|
}
|
|
|
|
Standard_Real mean = 0.5*theValue;
|
|
|
|
Standard_Real asc = Abs(fc-mean)*theKronrodW.Value(aNPnt2);
|
|
for(i = 1; i < aNPnt2; ++i) {
|
|
asc += theKronrodW.Value(i)*(Abs(f1(i) - mean) + Abs(f2(i) - mean));
|
|
}
|
|
|
|
asc *= aXr;
|
|
theValue *= aXr;
|
|
aGaussVal *= aXr;
|
|
|
|
// Compute the error and the new number of Kronrod points.
|
|
|
|
theError = Abs(theValue - aGaussVal);
|
|
|
|
Standard_Real scale =1.;
|
|
if(asc != 0. && theError != 0.) scale = Pow((200.*theError/asc), 1.5);
|
|
if(scale < 1.) theError = Min(theError, asc*scale);
|
|
|
|
//theFunction.GetStateNumber();
|
|
|
|
return Standard_True;
|
|
|
|
}
|
|
|