diff --git a/applications/utilities/miscellaneous/findRoot/findRoot.C b/applications/utilities/miscellaneous/findRoot/findRoot.C index 376c64f74..a6fe4d328 100644 --- a/applications/utilities/miscellaneous/findRoot/findRoot.C +++ b/applications/utilities/miscellaneous/findRoot/findRoot.C @@ -34,6 +34,7 @@ Description #include "IOmanip.H" #include "BisectionRoot.H" #include "RiddersRoot.H" +#include "NewtonSecantRoot.H" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // @@ -50,15 +51,38 @@ public: } }; + +class testFunctionDerivative +{ +public: + + testFunctionDerivative() + {} + + scalar operator()(const scalar& x) const + { + return 0.5*Foam::pow(x,-1.5) - Foam::cos(x); + } +}; + + int main(int argc, char *argv[]) { - testFunction tf; + testFunction tf; + testFunctionDerivative df; Info<< setprecision(10) << "Bisection root " << BisectionRoot(tf, 1e-5).root(0, 10) << nl << "Ridders root " - << RiddersRoot(tf, 1e-5).root(0, 10) << endl; + << RiddersRoot(tf, 1e-5).root(0, 10) << nl + << "NewtonSecant root " + << NewtonSecantRoot + ( + tf, + df, + 1e-5 + ).root(1.5) << endl; Info<< "End\n" << endl; diff --git a/src/ODE/findRoot/NewtonSecantRoot/NewtonSecantRoot.C b/src/ODE/findRoot/NewtonSecantRoot/NewtonSecantRoot.C new file mode 100644 index 000000000..96b5ea82e --- /dev/null +++ b/src/ODE/findRoot/NewtonSecantRoot/NewtonSecantRoot.C @@ -0,0 +1,105 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright held by original author + \\/ M anipulation | +------------------------------------------------------------------------------- +License + This file is part of OpenFOAM. + + OpenFOAM is free software; you can redistribute it and/or modify it + under the terms of the GNU General Public License as published by the + Free Software Foundation; either version 2 of the License, or (at your + option) any later version. + + OpenFOAM is distributed in the hope that it will be useful, but WITHOUT + ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License + for more details. + + You should have received a copy of the GNU General Public License + along with OpenFOAM; if not, write to the Free Software Foundation, + Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA + +Class + NewtonSecantRoot + +Description + NewtonSecant root finder for better stability. + Function is provided as a template parameter function object, evaluated + using operator()(const scalar x) + +Author + Aleksandar Jemcov. All rights reserved. + +\*---------------------------------------------------------------------------*/ + +#include "NewtonSecantRoot.H" + +// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // + +template +const Foam::label Foam::NewtonSecantRoot::maxIter = 60; + + +// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // + +template +Foam::NewtonSecantRoot::NewtonSecantRoot +( + const Func& f, + const Deriv& d, + const scalar eps +) +: + f_(f), + d_(d), + eps_(eps) +{} + + +// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // + +template +Foam::scalar Foam::NewtonSecantRoot::root +( + scalar xN +) const +{ + if (0 == d_(xN)) + { + FatalErrorIn + ( + "Foam::scalar Foam::NewtonSecantRoot::root\n" + "(\n" + " const scalar xN,\n" + ") const" + ) << "Jacobian equal to zero. f'(xN) = " << d_(xN) + << abort(FatalError); + } + + scalar xNp1; + + for (label nIter = 0; nIter < maxIter; ++nIter) + { + scalar fN = this->f_(xN); + scalar dN = this->d_(xN); + scalar dx = fN/dN; + + scalar xBar = xN - dx; + + xNp1 = xN - fN*fN/(dN*(fN - f_(xBar))); + if (mag(xN - xNp1) <= this->eps_) + { + return xNp1; + } + + xN = xNp1; + } + + return xNp1; +} + + +// ************************************************************************* // diff --git a/src/ODE/findRoot/NewtonSecantRoot/NewtonSecantRoot.H b/src/ODE/findRoot/NewtonSecantRoot/NewtonSecantRoot.H new file mode 100644 index 000000000..8673a9e5e --- /dev/null +++ b/src/ODE/findRoot/NewtonSecantRoot/NewtonSecantRoot.H @@ -0,0 +1,108 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright held by original author + \\/ M anipulation | +------------------------------------------------------------------------------- +License + This file is part of OpenFOAM. + + OpenFOAM is free software; you can redistribute it and/or modify it + under the terms of the GNU General Public License as published by the + Free Software Foundation; either version 2 of the License, or (at your + option) any later version. + + OpenFOAM is distributed in the hope that it will be useful, but WITHOUT + ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License + for more details. + + You should have received a copy of the GNU General Public License + along with OpenFOAM; if not, write to the Free Software Foundation, + Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA + +Class + NewtonSecantRoot + +Description + NewtonSecant root is based on combination of Newton and Secant methods. + Function and the funciton derivative are provided as a template parameter + function object, evaluated using operator()(const scalar x) + +Author + Aleksandar Jemcov. All rights reserved. + +SourceFiles + NewtonSecantRoot.C + +\*---------------------------------------------------------------------------*/ + +#ifndef NewtonSecantRoot_H +#define NewtonSecantRoot_H + +#include "scalar.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +/*---------------------------------------------------------------------------*\ + Class NewtonSecantRoot Declaration +\*---------------------------------------------------------------------------*/ + +template +class NewtonSecantRoot +{ + // Private data + + //- Reference to a function + const Func& f_; + + //- Reference to a derivative + const Deriv& d_; + + //- Required accuracy + const scalar eps_; + + +public: + + // Static data members + + //- Maximum number of iterations + static const label maxIter; + + + // Constructors + + //- Construct given a function + NewtonSecantRoot(const Func& f, const Deriv& d, const scalar eps); + + + // Destructor - default + + + // Member Functions + + //- Return root + scalar root(scalar xN) const; +}; + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#ifdef NoRepository +# include "NewtonSecantRoot.C" +#endif + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* //