Theoretica
A C++ numerical and automatic mathematical library
Loading...
Searching...
No Matches
roots.h
Go to the documentation of this file.
1
5
6#ifndef THEORETICA_ROOTS_H
7#define THEORETICA_ROOTS_H
8
9#include "../core/function.h"
10#include "../calculus/deriv.h"
11#include "../autodiff/dual.h"
12#include "../autodiff/dual2.h"
13#include "../algebra/vec.h"
14#include "../complex/complex.h"
15
16
17namespace theoretica {
18
19
27 template<typename RealFunction, typename Vector = vec2>
28 inline std::vector<Vector> find_root_intervals(
29 RealFunction f, real a, real b, unsigned int steps = 10) {
30
31 std::vector<Vector> res;
32 const real dx = (b - a) / (real) steps;
33
34 for (unsigned int i = 0; i < steps; ++i) {
35
36 const real x1 = a + i * dx;
37 const real x2 = a + (i + 1) * dx;
38
39 if(f(x1) * f(x2) <= 0)
40 res.push_back({x1, x2});
41 }
42
43 return res;
44 }
45
46
57 template<typename RealFunction>
59 RealFunction f, real a, real b,
60 real tol = OPTIMIZATION_TOL,
61 unsigned int max_iter = OPTIMIZATION_BISECTION_ITER) {
62
63 if(a > b) {
64 TH_MATH_ERROR("root_bisect", a, INVALID_ARGUMENT);
65 return nan();
66 }
67
68 if(f(a) * f(b) >= 0) {
69 TH_MATH_ERROR("root_bisect", f(a) * f(b), INVALID_ARGUMENT);
70 return nan();
71 }
72
73 real x_avg = 0.0;
74 real x_min = a;
75 real x_max = b;
76
77 unsigned int iter = 0;
78
79 while((x_max - x_min) > (2 * tol) && iter <= max_iter) {
80
81 x_avg = (x_max + x_min) / 2.0;
82
83 if(f(x_avg) * f(x_min) > 0)
84 x_min = x_avg;
85 else
86 x_max = x_avg;
87
88 ++iter;
89 }
90
91 if(iter > max_iter) {
92 TH_MATH_ERROR("root_bisect", x_avg, NO_ALGO_CONVERGENCE);
93 return nan();
94 }
95
96 return x_avg;
97 }
98
99
119 template<typename RealFunction>
121 RealFunction f, real a, real b, real tol = OPTIMIZATION_TOL,
122 unsigned int n0 = 1, real k1 = 0.0) {
123
124 // Default value for k1
125 if (k1 == 0.0)
126 k1 = 0.2 / (b - a);
127
128 if(a > b) {
129 TH_MATH_ERROR("root_itp", a, INVALID_ARGUMENT);
130 return nan();
131 }
132
133 real y_a = f(a);
134 real y_b = f(b);
135
136 if(y_a * y_b >= 0) {
137 TH_MATH_ERROR("root_itp", y_a * y_b, INVALID_ARGUMENT);
138 return nan();
139 }
140
141 // Monotonicity of the function
142 const int monotone = (y_a < y_b) ? 1 : -1;
143
144 real x_t, x_new;
145
146 const long int n_half = th::floor(th::log2((b - a) / tol));
147 const long int n_max = n_half + n0;
148
149 real eps = tol * pow(2.0, n_max);
150 long int iter = 0;
151
152 while((b - a) > (2 * tol) && iter <= n_max) {
153
154
155 // Interpolation
156 const real x_f = (a * y_b - b * y_a) / (y_b - y_a);
157 const real x_half = 0.5 * (a + b);
158
159
160 // Truncation
161 const int sigma = sgn(x_half - x_f);
162 const real delta = k1 * square(b - a);
163
164 if (delta <= abs(x_half - x_f))
165 x_t = x_f + sigma * delta;
166 else
167 x_t = x_half;
168
169
170 // Projection
171 const real r = eps - (b - a) / 2.0;
172
173 if (abs(x_t - x_half) <= r)
174 x_new = x_t;
175 else
176 x_new = x_half - sigma * r;
177
178
179 // Update
180 const real y_new = f(x_new);
181
182 if (monotone * y_new > 0) {
183 b = x_new;
184 y_b = y_new;
185 } else if (monotone * y_new < 0) {
186 a = x_new;
187 y_a = y_new;
188 } else {
189 return x_new;
190 }
191
192 eps *= 0.5;
193 iter++;
194 }
195
196 if(abs(b - a) > 2 * tol) {
197 TH_MATH_ERROR("root_itp", (a + b) / 2.0, NO_ALGO_CONVERGENCE);
198 return nan();
199 }
200
201 return (a + b) / 2.0;
202 }
203
204
216 template<typename RealFunction>
219 real tol = OPTIMIZATION_TOL, unsigned int max_iter = OPTIMIZATION_NEWTON_ITER) {
220
221 real x = guess;
222 real f_x = inf();
223 unsigned int iter = 0;
224
225 while(abs(f_x) > tol && iter <= max_iter) {
226
227 f_x = f(x);
228 x = x - (f_x / Df(x));
229 iter++;
230 }
231
232 if(iter > max_iter) {
233 TH_MATH_ERROR("root_newton", x, NO_ALGO_CONVERGENCE);
234 return nan();
235 }
236
237 return x;
238 }
239
240
254 dual(*f)(dual), real guess = 0,
255 real tol = OPTIMIZATION_TOL,
256 unsigned int max_iter = OPTIMIZATION_NEWTON_ITER) {
257
258 real x = guess;
259 unsigned int iter = 0;
260
261 dual s = dual(inf(), 0.0);
262
263 while(abs(s.Re()) > tol && iter <= max_iter) {
264
265 // Compute the function and its derivative at the same time
266 s = f(dual(x, 1.0));
267
268 x = x - (s.Re() / s.Dual());
269 iter++;
270 }
271
272 if(iter > max_iter) {
273 TH_MATH_ERROR("root_newton", x, NO_ALGO_CONVERGENCE);
274 return nan();
275 }
276
277 return x;
278 }
279
280
292 template <
293 typename Type = real,
294 typename ComplexFunction = std::function<complex<Type>(complex<Type>)>
295 >
298 real tol = OPTIMIZATION_TOL,
299 unsigned int max_iter = OPTIMIZATION_NEWTON_ITER) {
300
303 unsigned int iter = 0;
304
305 while(f_z * f_z.conjugate() > tol * tol && iter <= max_iter) {
306
307 f_z = f(z);
308 z = z - (f_z / df(z));
309 iter++;
310 }
311
312 if(iter > max_iter) {
313 TH_MATH_ERROR("root_newton", z.Re(), NO_ALGO_CONVERGENCE);
314 return complex<Type>(nan(), nan());
315 }
316
317 return z;
318 }
319
320
333 template<typename RealFunction>
337 real tol = OPTIMIZATION_TOL,
338 unsigned int max_iter = OPTIMIZATION_HALLEY_ITER) {
339
340 real x = guess;
341 real f_x = inf();
342 unsigned int iter = 0;
343
344 while(abs(f_x) > tol && iter <= max_iter) {
345
346 f_x = f(x);
347 const real df_x = Df(x);
348
349 x = x - (2 * f_x * df_x) / (2 * square(df_x) - f_x * D2f(x));
350 iter++;
351 }
352
353 if(iter > max_iter) {
354 TH_MATH_ERROR("root_halley", x, NO_ALGO_CONVERGENCE);
355 return nan();
356 }
357
358 return x;
359 }
360
361
376 dual2(*f)(dual2), real guess = 0,
377 real tol = OPTIMIZATION_TOL,
378 unsigned int max_iter = OPTIMIZATION_HALLEY_ITER) {
379
380 real x = guess;
381 dual2 s = dual2(inf(), 0.0, 0.0);
382 unsigned int iter = 0;
383
384 while(abs(s.Re()) > tol && iter <= max_iter) {
385
386 // Compute the function value and the first and
387 // second derivatives at the same time.
388 s = f(dual2(x, 1, 0));
389
390 const real f_x = s.Re();
391 const real df_x = s.Dual1();
392 const real d2f_x = s.Dual2();
393
394 x = x - (2 * f_x * df_x) / (2 * square(df_x) - f_x * d2f_x);
395 iter++;
396 }
397
398 if(iter > max_iter) {
399 TH_MATH_ERROR("root_halley", x, NO_ALGO_CONVERGENCE);
400 return nan();
401 }
402
403 return x;
404 }
405
406
417 template<typename RealFunction>
419 RealFunction f, real guess = 0,
420 real tol = OPTIMIZATION_TOL,
421 unsigned int max_iter = OPTIMIZATION_STEFFENSEN_ITER) {
422
423 real x = guess;
424 real f_x = inf();
425 unsigned int iter = 0;
426
427 while(abs(f_x) > tol && iter <= max_iter) {
428
429 const real f_x = f(x);
430 const real g_x = (f(x + f_x) / f_x) - 1;
431
432 x = x - (f_x / g_x);
433 iter++;
434 }
435
436 if(iter > max_iter) {
437 TH_MATH_ERROR("root_steffensen", x, NO_ALGO_CONVERGENCE);
438 return nan();
439 }
440
441 return x;
442 }
443
444
462 template<typename RealFunction>
466 real tol = OPTIMIZATION_TOL,
467 unsigned int max_iter = OPTIMIZATION_CHEBYSHEV_ITER) {
468
469 real x = guess;
470 real f_x = inf();
471 unsigned int iter = 0;
472
473 while(abs(f_x) > tol && iter <= max_iter) {
474
475 f_x = f(x);
476 const real df_x = Df(x);
477 const real u = f_x / df_x;
478
479 x = x - u - square(u) * D2f(x) / (2.0 * df_x);
480
481 iter++;
482 }
483
484 if(iter > max_iter) {
485 TH_MATH_ERROR("root_chebyshev", x, NO_ALGO_CONVERGENCE);
486 return nan();
487 }
488
489 return x;
490 }
491
492
512 dual2(*f)(dual2), real guess = 0,
513 real tol = OPTIMIZATION_TOL,
514 unsigned int max_iter = OPTIMIZATION_CHEBYSHEV_ITER) {
515
516 real x = guess;
517 dual2 s = dual2(inf(), 0.0, 0.0);
518 unsigned int iter = 0;
519
520 while(abs(s.Re()) > tol && iter <= max_iter) {
521
522 s = f(dual2(x, 1.0, 0.0));
523
524 const real f_x = s.Re();
525 const real df_x = s.Dual1();
526 const real u = f_x / df_x;
527
528 x = x - u - square(u) * s.Dual2() / (2.0 * df_x);
529 iter++;
530 }
531
532 if(iter > max_iter) {
533 TH_MATH_ERROR("root_chebyshev", x, NO_ALGO_CONVERGENCE);
534 return nan();
535 }
536
537 return x;
538 }
539
540
557 template<typename RealFunction>
560 real tol = OPTIMIZATION_TOL,
561 unsigned int max_iter = OPTIMIZATION_OSTROWSKI_ITER) {
562
563 real x = guess;
564 real f_x = inf();
565 unsigned int iter = 0;
566
567 while(abs(f_x) > tol && iter <= max_iter) {
568
569 f_x = f(x);
570 const real df_x = Df(x);
571 const real u = f_x / df_x;
572 const real f_xu = f(x - u);
573
574 x = x - u - (f_xu / df_x) * (f_x / (f_x - 2 * f_xu));
575
576 iter++;
577 }
578
579 if(iter > max_iter) {
580 TH_MATH_ERROR("root_ostrowski", x, NO_ALGO_CONVERGENCE);
581 return nan();
582 }
583
584 return x;
585 }
586
587
605 template<typename RealFunction>
608 real tol = OPTIMIZATION_TOL,
609 unsigned int max_iter = OPTIMIZATION_JARRAT_ITER) {
610
611 real x = guess;
612 real f_x = inf();
613 unsigned int iter = 0;
614
615 while(abs(f_x) > tol && iter <= max_iter) {
616
617 f_x = f(x);
618 const real df_x = Df(x);
619 const real u = f_x / df_x;
620 const real f_xu = Df(x - 2.0 * u / 3.0);
621
622 x = x - 0.5 * u + f_x / (df_x - 3 * f_xu);
623 iter++;
624 }
625
626 if(iter > max_iter) {
627 TH_MATH_ERROR("root_jarrat", x, NO_ALGO_CONVERGENCE);
628 return nan();
629 }
630
631 return x;
632 }
633
634
649 template<typename RealFunction>
650 inline std::vector<real> roots(
651 RealFunction f, real a, real b,
652 real tol = OPTIMIZATION_TOL, real steps = 10) {
653
654 if(steps == 0) {
655 TH_MATH_ERROR("roots", steps, DIV_BY_ZERO);
656 return {nan()};
657 }
658
659 // Find candidate intervals
660 std::vector<vec2> intervals = find_root_intervals(f, a, b, steps);
661
662 std::vector<real> res;
663 res.reserve(intervals.size());
664
665 // Iterate through all candidate intervals and refine the results
666 for (unsigned int i = 0; i < intervals.size(); ++i) {
667
668 // Check whether the extremes of the candidate intervals
669 // happen to coincide with the roots
670 if(abs(f(intervals[i][0])) <= MACH_EPSILON) {
671 res.push_back(intervals[i][0]);
672 continue;
673 }
674
675 if(abs(f(intervals[i][1])) <= MACH_EPSILON) {
676 res.push_back(intervals[i][1]);
677 continue;
678 }
679
680 // Approximate the roots using bisection inside each interval
681 res.push_back(
682 root_bisect(f, intervals[i][0], intervals[i][1], tol)
683 );
684 }
685
686 return res;
687 }
688
689
697 template<typename Field>
698 inline std::vector<Field> roots(
699 const polynomial<Field>& p,
700 real tolerance = OPTIMIZATION_TOL,
701 unsigned int steps = 0) {
702
703 // Effective order of the polynomial
704 const unsigned int n = p.find_order();
705 p /= p.coeff[n];
706
707 // Absolute value of the highest coefficient
708 Field a_hi = abs(p.coeff[n]);
709 Field a_sum = 0;
710
711 // Sum the absolute values of the lesser coefficients
712 for (unsigned int i = 0; i < n; ++i)
713 a_sum += abs(p.coeff[i]);
714
715 // The roots are bounded in absolute value by the maximum
716 const Field M = max(a_hi, a_sum);
717
718 // Back track from the bounds to the first sign inversion ?
719
720 return roots(
721 [p](real x) { return p(x); },
722 -M, M, tolerance, steps ? steps : (2 * n));
723 }
724
725}
726
727#endif
Second order dual number class.
Definition dual2.h:29
Dual number class.
Definition dual.h:28
#define TH_MATH_ERROR(F_NAME, VALUE, EXCEPTION)
TH_MATH_ERROR is a macro which throws exceptions or modifies errno (depending on which compiling opti...
Definition error.h:219
Main namespace of the library which contains all functions and objects.
Definition algebra.h:27
double real
A real number, defined as a floating point type.
Definition constants.h:198
real root_newton(RealFunction f, RealFunction Df, real guess=0.0, real tol=OPTIMIZATION_TOL, unsigned int max_iter=OPTIMIZATION_NEWTON_ITER)
Find a root of a univariate real function using Newton's method.
Definition roots.h:217
real root_halley(RealFunction f, RealFunction Df, RealFunction D2f, real guess=0, real tol=OPTIMIZATION_TOL, unsigned int max_iter=OPTIMIZATION_HALLEY_ITER)
Find a root of a univariate real function using Halley's method.
Definition roots.h:334
real inf()
Return positive infinity in floating point representation.
Definition error.h:76
dual2 log2(dual2 x)
Compute the natural logarithm of a second order dual number.
Definition dual2_functions.h:166
std::remove_reference_t< decltype(std::declval< Structure >()[0])> vector_element_t
Extract the type of a vector (or any indexable container) from its operator[].
Definition core_traits.h:134
std::vector< Vector > find_root_intervals(RealFunction f, real a, real b, unsigned int steps=10)
Find candidate intervals for root finding by evaluating a function at equidistant points inside an in...
Definition roots.h:28
std::vector< real > roots(RealFunction f, real a, real b, real tol=OPTIMIZATION_TOL, real steps=10)
Find the roots of a univariate real function inside a given interval, by first searching for candidat...
Definition roots.h:650
auto max(const Vector &X)
Finds the maximum value inside a dataset.
Definition dataset.h:330
real root_bisect(RealFunction f, real a, real b, real tol=OPTIMIZATION_TOL, unsigned int max_iter=OPTIMIZATION_BISECTION_ITER)
Find the root of a univariate real function using bisection inside a compact interval [a,...
Definition roots.h:58
constexpr real MACH_EPSILON
Machine epsilon for the real type.
Definition constants.h:207
real root_itp(RealFunction f, real a, real b, real tol=OPTIMIZATION_TOL, unsigned int n0=1, real k1=0.0)
Find a root of a univariate real function using the ITP (Interpolate-Truncate-Project) method,...
Definition roots.h:120
real root_steffensen(RealFunction f, real guess=0, real tol=OPTIMIZATION_TOL, unsigned int max_iter=OPTIMIZATION_STEFFENSEN_ITER)
Find a root of a univariate real function using Steffensen's method.
Definition roots.h:418
real root_ostrowski(RealFunction f, RealFunction Df, real guess=0.0, real tol=OPTIMIZATION_TOL, unsigned int max_iter=OPTIMIZATION_OSTROWSKI_ITER)
Find a root of a univariate real function using Ostrowski's method.
Definition roots.h:558
real root_chebyshev(RealFunction f, RealFunction Df, RealFunction D2f, real guess=0, real tol=OPTIMIZATION_TOL, unsigned int max_iter=OPTIMIZATION_CHEBYSHEV_ITER)
Find a root of a univariate real function using Chebyshev's method.
Definition roots.h:463
int sgn(real x)
Return the sign of x (1 if positive, -1 if negative, 0 if null)
Definition real_analysis.h:259
real nan()
Return a quiet NaN number in floating point representation.
Definition error.h:54
real root_jarrat(RealFunction f, RealFunction Df, real guess=0.0, real tol=OPTIMIZATION_TOL, unsigned int max_iter=OPTIMIZATION_JARRAT_ITER)
Find a root of a univariate real function using Jarrat's method.
Definition roots.h:606
TH_CONSTEXPR int floor(real x)
Compute the floor of x Computes the maximum integer number that is smaller than x.
Definition real_analysis.h:271