include/opencascade/MathRoot_Secant.hxx

0001 // Copyright (c) 2025 OPEN CASCADE SAS
0002 //
0003 // This file is part of Open CASCADE Technology software library.
0004 //
0005 // This library is free software; you can redistribute it and/or modify it under
0006 // the terms of the GNU Lesser General Public License version 2.1 as published
0007 // by the Free Software Foundation, with special exception defined in the file
0008 // OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT
0009 // distribution for complete text of the license and disclaimer of any warranty.
0010 //
0011 // Alternatively, this file may be used under the terms of Open CASCADE
0012 // commercial license or contractual agreement.
0013
0014 #ifndef _MathRoot_Secant_HeaderFile
0015 #define _MathRoot_Secant_HeaderFile
0016
0017 #include <MathUtils_Types.hxx>
0018 #include <MathUtils_Config.hxx>
0019 #include <MathUtils_Core.hxx>
0020 #include <MathUtils_Convergence.hxx>
0021
0022 #include <cmath>
0023
0024 namespace MathRoot
0025 {
0026 using namespace MathUtils;
0027
0028 //! Secant method for root finding.
0029 //! Does not require derivative, uses finite difference approximation.
0030 //! Converges superlinearly (order ~1.618, the golden ratio).
0031 //!
0032 //! Algorithm:
0033 //! 1. Start with two initial points x0, x1
0034 //! 2. Approximate derivative: f'(x) ~ (f(x1) - f(x0)) / (x1 - x0)
0035 //! 3. Newton-like update: x2 = x1 - f(x1) * (x1 - x0) / (f(x1) - f(x0))
0036 //! 4. Repeat with x0 = x1, x1 = x2
0037 //!
0038 //! @tparam Function type with Value(double theX, double& theF) method
0039 //! @param theFunc function object providing only value
0040 //! @param theX0 first initial point
0041 //! @param theX1 second initial point (different from theX0)
0042 //! @param theConfig solver configuration
0043 //! @return result containing root location and convergence status
0044 template <typename Function>
0045 MathUtils::ScalarResult Secant(Function&                theFunc,
0046                                double                   theX0,
0047                                double                   theX1,
0048                                const MathUtils::Config& theConfig = MathUtils::Config())
0049 {
0050   MathUtils::ScalarResult aResult;
0051
0052   double aX0 = theX0;
0053   double aX1 = theX1;
0054   double aF0 = 0.0;
0055   double aF1 = 0.0;
0056
0057   // Evaluate at initial points
0058   if (!theFunc.Value(aX0, aF0))
0059   {
0060     aResult.Status = MathUtils::Status::NumericalError;
0061     return aResult;
0062   }
0063   if (!theFunc.Value(aX1, aF1))
0064   {
0065     aResult.Status = MathUtils::Status::NumericalError;
0066     return aResult;
0067   }
0068
0069   for (int anIter = 0; anIter < theConfig.MaxIterations; ++anIter)
0070   {
0071     aResult.NbIterations = anIter + 1;
0072
0073     // Check convergence
0074     if (MathUtils::IsFConverged(aF1, theConfig.FTolerance))
0075     {
0076       aResult.Status = MathUtils::Status::OK;
0077       aResult.Root   = aX1;
0078       aResult.Value  = aF1;
0079       return aResult;
0080     }
0081
0082     // Secant step
0083     const double aDenom = aF1 - aF0;
0084     if (MathUtils::IsZero(aDenom))
0085     {
0086       aResult.Status = MathUtils::Status::NumericalError;
0087       aResult.Root   = aX1;
0088       aResult.Value  = aF1;
0089       return aResult;
0090     }
0091
0092     const double aXNew = aX1 - aF1 * (aX1 - aX0) / aDenom;
0093
0094     // Check X convergence
0095     if (MathUtils::IsXConverged(aX1, aXNew, theConfig.XTolerance))
0096     {
0097       double aFNew = 0.0;
0098       theFunc.Value(aXNew, aFNew);
0099       aResult.Status = MathUtils::Status::OK;
0100       aResult.Root   = aXNew;
0101       aResult.Value  = aFNew;
0102       return aResult;
0103     }
0104
0105     // Update for next iteration
0106     aX0 = aX1;
0107     aF0 = aF1;
0108     aX1 = aXNew;
0109
0110     if (!theFunc.Value(aX1, aF1))
0111     {
0112       aResult.Status = MathUtils::Status::NumericalError;
0113       aResult.Root   = aX1;
0114       return aResult;
0115     }
0116   }
0117
0118   aResult.Status = MathUtils::Status::MaxIterations;
0119   aResult.Root   = aX1;
0120   aResult.Value  = aF1;
0121   return aResult;
0122 }
0123
0124 //! Secant method with automatic initial points.
0125 //! Creates second point by small perturbation of the initial guess.
0126 //!
0127 //! @tparam Function type with Value(double theX, double& theF) method
0128 //! @param theFunc function object providing only value
0129 //! @param theX0 initial guess
0130 //! @param theConfig solver configuration
0131 //! @return result containing root location and convergence status
0132 template <typename Function>
0133 MathUtils::ScalarResult SecantAuto(Function&                theFunc,
0134                                    double                   theX0,
0135                                    const MathUtils::Config& theConfig = MathUtils::Config())
0136 {
0137   // Create second point by small perturbation
0138   const double aDelta = (std::abs(theX0) > 1.0) ? 0.01 * theX0 : 0.01;
0139   return Secant(theFunc, theX0, theX0 + aDelta, theConfig);
0140 }
0141
0142 } // namespace MathRoot
0143
0144 #endif // _MathRoot_Secant_HeaderFile