suanPan
ridders.hpp
Go to the documentation of this file.
1/*******************************************************************************
2 * Copyright (C) 2017-2024 Theodore Chang
3 *
4 * This program is free software: you can redistribute it and/or modify
5 * it under the terms of the GNU General Public License as published by
6 * the Free Software Foundation, either version 3 of the License, or
7 * (at your option) any later version.
8 *
9 * This program is distributed in the hope that it will be useful,
10 * but WITHOUT ANY WARRANTY; without even the implied warranty of
11 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12 * GNU General Public License for more details.
13 *
14 * You should have received a copy of the GNU General Public License
15 * along with this program. If not, see <http://www.gnu.org/licenses/>.
16 ******************************************************************************/
17
18#ifndef RIDDERS_HPP
19#define RIDDERS_HPP
20
21#include <suanPan.h>
22
23template<std::invocable<double> T> double ridders(const T& func, double x1, double f1, double x2, double f2, const double tolerance) {
24 double target;
25
26 auto counter = 0u;
27 while(true) {
28 counter += 2u;
29
30 const auto x3 = .5 * (x1 + x2);
31 const auto f3 = func(x3);
32 if(std::fabs(f3) < tolerance || fabs(x2 - x1) < tolerance) {
33 target = x3;
34 break;
35 }
36
37 const auto dx = (x3 - x1) * f3 / std::sqrt(f3 * f3 - f1 * f2);
38
39 const auto x4 = f1 > f2 ? x3 + dx : x3 - dx;
40 const auto f4 = func(x4);
41 if(std::fabs(f4) < tolerance) {
42 target = x4;
43 break;
44 }
45
46 // one end is x4
47 // pick the other from x3, x2, x1
48 if(f4 * f3 < 0.) {
49 x1 = x3;
50 f1 = f3;
51 }
52 else if(f4 * f2 < 0.) {
53 x1 = x2;
54 f1 = f2;
55 }
56
57 x2 = x4;
58 f2 = f4;
59 }
60
61 suanpan_debug("Ridders' method initial guess {:.5E} with {} iterations.\n", target, counter);
62
63 return target;
64}
65
66#endif
67
double ridders(const T &func, double x1, double f1, double x2, double f2, const double tolerance)
Definition: ridders.hpp:23
#define suanpan_debug(...)
Definition: suanPan.h:307