Secant method in Python
Step to the zero of the line through the last two points. Order (1 + 5^(1/2))/2 ≈ 1.618 near a simple root, with one evaluation of f per step and no derivative.
- File
secant_method.py, 62 lines- Level
- I0 not checked yet
- Toolchain
- Python 3.14, pytest 9
- Test framework
- pytest
- Run the tests
python3 -m pytest
# SPDX-License-Identifier: MIT
# Part of https://algorithm.tamnd.com/secant-method/. The authors are listed in prov/8VVH8M.jsonl.
"""The secant method: step to the zero of the line through the last two points."""
import cmath
import math
from dataclasses import dataclass
from typing import Any, Callable
@dataclass
class SecantResult:
"""The returned point, why the loop ended, and the work done."""
root: Any
reason: str # "root", "converged", "failure" or "kmax reached"
iterations: int
evaluations: int
@property
def converged(self) -> bool:
return self.reason in ("root", "converged")
def _finite(v: Any) -> bool:
if isinstance(v, complex):
return cmath.isfinite(v)
return math.isfinite(float(v))
def secant(f: Callable[[Any], Any], x0: Any, x1: Any, tol_abs: float = 0.0,
tol_rel: float = 2.0 ** -51, max_iter: int = 100) -> SecantResult:
"""Run the secant method on f from the two points x0 != x1.
Works for float, complex, Decimal and Fraction values. Precondition:
tol_abs >= 0, tol_rel >= 0, not both zero. Convergence is guaranteed only
near a simple root of a C^2 function. Cost: at most max_iter + 2 calls of
f, O(1) extra space.
"""
f0 = f(x0) # @line L1
if f0 == 0:
return SecantResult(x0, "root", 0, 1)
f1 = f(x1) # @line L2
evals = 2
if f1 == 0:
return SecantResult(x1, "root", 0, evals)
for k in range(1, max_iter + 1): # @line L3
if f1 == f0 or not _finite(f0) or not _finite(f1): # @line L4
return SecantResult(x1, "failure", k, evals)
s = -f1 * (x1 - x0) / (f1 - f0) # @line L5
x2 = x1 + s # @line L6
if not _finite(x2): # @line L7
return SecantResult(x1, "failure", k, evals)
if abs(s) <= tol_abs + tol_rel * abs(x2): # @line L8
return SecantResult(x2, "converged", k, evals)
f2 = f(x2) # @line L9
evals += 1
if f2 == 0:
return SecantResult(x2, "root", k, evals)
x0, f0, x1, f1 = x1, f1, x2, f2 # @line L10
return SecantResult(x1, "kmax reached", max_iter, evals) # @line L11
test_secant_method.py 110 lines
# SPDX-License-Identifier: MIT
# Part of https://algorithm.tamnd.com/secant-method/. The authors are listed in prov/8VVH8M.jsonl.
import math
import random
from decimal import Decimal, getcontext
from secant_method import secant
PHI = (1 + 5 ** 0.5) / 2
def test_worked_example_table():
xs = []
def f(x):
xs.append(x)
return math.cos(x) - x
r = secant(f, 0.0, 1.0)
assert xs == [0.0, 1.0, 0.6850733573260451, 0.736298997613654, 0.7391193619116293,
0.7390851121274639, 0.7390851332150012, 0.7390851332151607]
assert r.reason == "root" and r.evaluations == 8 and r.root == 0.7390851332151607
def test_known_roots():
cases = [
(lambda x: math.exp(x) - 2, 0.0, 1.0, math.log(2)),
(lambda x: x ** 3 - x - 2, 1.0, 2.0, 1.5213797068045676),
(lambda x: (x - 1) * (x + 2) * (x - 3.5), 3.0, 4.0, 3.5),
]
for f, a, b, root in cases:
r = secant(f, a, b)
assert r.converged and abs(r.root - root) <= 4e-16 * max(1.0, abs(root))
def test_complex_root():
r = secant(lambda z: z * z + 1, complex(0.5, 0.5), complex(0.6, 0.8))
assert r.converged and abs(r.root - 1j) < 1e-14
def test_observed_order_is_golden_ratio():
getcontext().prec = 400
root = Decimal("1.5213797068045675696") # refined below
for _ in range(12):
root -= (root ** 3 - root - 2) / (3 * root * root - 1)
xs = []
def f(x):
xs.append(x)
return x ** 3 - x - 2
secant(f, Decimal(1), Decimal(2), tol_abs=Decimal("1e-300"), tol_rel=Decimal(0),
max_iter=11)
errs = [abs(x - root) for x in xs]
qs = [float((errs[k + 1] / errs[k]).ln() / (errs[k] / errs[k - 1]).ln())
for k in range(7, len(errs) - 1)]
assert all(abs(q - PHI) < 0.01 for q in qs)
c = 6 * root / (2 * (3 * root * root - 1))
assert abs(float(errs[-1] / errs[-2] ** Decimal(PHI)) - float(c) ** (1 / PHI)) < 1e-2
def test_double_root_is_linear():
getcontext().prec = 60
xs = []
def f(x):
xs.append(x)
return x ** 3 - 3 * x + 2
secant(f, Decimal(2), Decimal("1.5"), tol_abs=Decimal("1e-15"), tol_rel=Decimal(0),
max_iter=60)
ratios = [float((xs[k + 1] - 1) / (xs[k] - 1)) for k in range(25, 30)]
assert all(abs(q - 1 / PHI) < 1e-3 for q in ratios)
def test_evaluation_counts():
r = secant(lambda x: x - 0.25, 0.0, 1.0)
# the line is exact: x2 = 0.25 is a root found at the first pass
assert r.reason == "root" and r.evaluations == 3 and r.iterations == 1
r = secant(lambda x: x ** 3 - x - 2, 1.0, 2.0, tol_rel=1e-6)
assert r.reason == "converged" and r.evaluations == r.iterations + 1
def test_edges():
assert secant(lambda x: x, 0.0, 1.0).reason == "root"
assert secant(lambda x: x - 1.0, 0.0, 1.0).reason == "root"
r = secant(lambda x: x * x - 1, -0.5, 0.5)
assert r.reason == "failure"
r = secant(lambda x: float("nan"), 0.0, 1.0)
assert r.reason == "failure"
r = secant(lambda x: math.atan(x), 2.0, 3.0)
assert not r.converged
r = secant(lambda x: x * x + 1, 0.0, 1.0, max_iter=5)
assert r.reason in ("kmax reached", "failure") and r.evaluations <= 7
def test_random_cubics_with_known_roots():
rng = random.Random(1618)
count = 0
while count < 500:
roots = sorted(rng.uniform(-5, 5) for _ in range(3))
if roots[1] - roots[0] < 1 or roots[2] - roots[1] < 1:
continue
count += 1
a, b, c = roots
target = rng.choice(roots)
x0 = target + rng.uniform(-0.05, 0.05)
r = secant(lambda x: (x - a) * (x - b) * (x - c), x0, x0 + 0.01)
assert r.converged and abs(r.root - target) <= 1e-13
What the next level needs
- the tests did not run in this build