Alexandria 2.31.0
SDC-CH common library for the Euclid project
Loading...
Searching...
No Matches
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
22namespace Euclid {
23namespace MathUtils {
24
25SecantReturn 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)) {
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;
54 break;
55 } else if (x2 >= params.max) {
56 ret.root = params.max;
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
Interface class representing a function with an arbitrary number of parameters.
Definition Function.h:104
T isfinite(T... args)
SecantReturn secantMethod(const Function &func, double x0, double x1, const SecantParams &params=SecantParams{})
std::size_t max_iter
Maximum number of iterations.
double min
If the gradient moves the next iteration below this limit, clip the result.
double max
If the gradient moves the next iteration above this limit, clip the result.