5.4.2 Root Finding (Iteration)

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85

/*

Finds an $x$ for a continuous function $f$ such that $f(x) = 0$ using iterative approximation by an
initial guess that is close to the answer. Each step follows the tangent line at the current guess
down to its $x$-intercept, that is, $x \leftarrow x - f(x)/f'(x)$. Newton's method requires an
explicit definition of the function's derivative while the secant method starts with two initial
guesses and approximates the derivative using the secant slope from the previous iteration. With
smooth functions and good initial guesses, Newton's method has quadratic convergence, while the
secant method has convergence order approximately $1.618$.

- `newton_root(f, fprime, x0, eps = 1e-15, iterations = 100)` returns a root $x$ for a function `f`
  with derivative `fprime` using an initial guess `x0` which should be relatively close to $x$.
- `secant_root(f, x0, x1, eps = 1e-15, iterations = 100)` returns a root $x$ for a function `f`
  using two initial guesses `x0` and `x1` which should be relatively close to $x$.

Time Complexity:
- O(n) calls to `f()` and `fprime()` per call to `newton_root()`, and O(n) calls to `f()` per call
  to `secant_root()`, where $n$ is the number of iterations.

Space Complexity:
- O(1) auxiliary for both operations.

*/

#include <cmath>
#include <stdexcept>

template<typename Fn, typename Deriv>
double newton_root(Fn f, Deriv fprime, double x0, double eps = 1e-15, int iterations = 100) {
  double x = x0, error = eps + 1;
  for (int i = 0; std::isfinite(error) && error > eps && i < iterations; i++) {
    double fx = f(x);
    if (fx == 0) {
      return x;
    }
    double xnew = x - fx / fprime(x);
    error = fabs(xnew - x);
    x = xnew;
  }
  if (!std::isfinite(error) || error > eps) {
    throw std::runtime_error("Newton's method failed to converge.");
  }
  return x;
}

template<typename Fn>
double secant_root(Fn f, double x0, double x1, double eps = 1e-15, int iterations = 100) {
  double xold = x0, fxold = f(x0), x = x1, error = eps + 1;
  for (int i = 0; std::isfinite(error) && error > eps && i < iterations; i++) {
    double fx = f(x);
    if (fx == 0) {
      return x;
    }
    double xnew = x - fx * ((x - xold) / (fx - fxold));
    error = fabs(xnew - x);
    xold = x;
    fxold = fx;
    x = xnew;
  }
  if (!std::isfinite(error) || error > eps) {
    throw std::runtime_error("Secant method failed to converge.");
  }
  return x;
}

/*** Example Usage ***/

#include <cassert>

double f(double x) {
  return x * x - 4 * sin(x);
}

double fprime(double x) {
  return 2 * x - 4 * cos(x);
}

int main() {
  assert(fabs(f(newton_root(f, fprime, 3))) < 1e-10);
  assert(fabs(f(secant_root(f, 3, 2))) < 1e-10);
  auto g = [](double x) { return x * x - 2; };
  auto gprime = [](double x) { return 2 * x; };
  assert(fabs(g(newton_root(g, gprime, 1))) < 1e-10);
  return 0;
}