Alexandria  2.27.0
SDC-CH common library for the Euclid project
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
SecantMethod.cpp
Go to the documentation of this file.
1 /*
2  * Copyright (C) 2022 Euclid Science Ground Segment
3  *
4  * This library is free software; you can redistribute it and/or modify it under
5  * the terms of the GNU Lesser General Public License as published by the Free
6  * Software Foundation; either version 3.0 of the License, or (at your option)
7  * any later version.
8  *
9  * This library is distributed in the hope that it will be useful, but WITHOUT
10  * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
11  * FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
12  * details.
13  *
14  * You should have received a copy of the GNU Lesser General Public License
15  * along with this library; if not, write to the Free Software Foundation, Inc.,
16  * 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
17  */
18 
20 #include <cmath>
21 
22 namespace Euclid {
23 namespace MathUtils {
24 
25 SecantReturn secantMethod(const Function& func, double x0, double x1, const SecantParams& params) {
27 
28  double y0 = func(x0);
29 
30  if (!std::isfinite(y0)) {
31  return {x0, SecantEndReason::VALUE_ERROR, 0};
32  }
33 
34  for (ret.iterations = 0; ret.iterations < params.max_iter; ++ret.iterations) {
35  if (std::abs(y0) <= params.atol) {
36  break;
37  }
38 
39  double y1 = func(x1);
40  if (!std::isfinite(y1)) {
41  ret.reason = SecantEndReason::VALUE_ERROR;
42  break;
43  }
44  if (y1 == y0) {
45  ret.reason = SecantEndReason::GRADIENT;
46  break;
47  }
48 
49  double x2 = x1 - y1 * ((x1 - x0) / (y1 - y0));
50 
51  if (x2 <= params.min) {
52  ret.root = params.min;
53  ret.reason = SecantEndReason::OUT_OF_BOUNDS;
54  break;
55  } else if (x2 >= params.max) {
56  ret.root = params.max;
57  ret.reason = SecantEndReason::OUT_OF_BOUNDS;
58  break;
59  }
60  x0 = x1;
61  x1 = x2;
62  y0 = y1;
63  ret.root = x0;
64  }
65  return ret;
66 }
67 
68 } // namespace MathUtils
69 } // namespace Euclid
SecantReturn secantMethod(const Function &func, double x0, double x1, const SecantParams &params=SecantParams{})
Interface class representing a function with an arbitrary number of parameters.
Definition: Function.h:104
double max
If the gradient moves the next iteration above this limit, clip the result.
Definition: SecantMethod.h:36
std::size_t max_iter
Maximum number of iterations.
Definition: SecantMethod.h:31
T isfinite(T...args)
double min
If the gradient moves the next iteration below this limit, clip the result.
Definition: SecantMethod.h:34