158 lines
5.0 KiB
Python
158 lines
5.0 KiB
Python
"""Stage 0 - does the similarity-proportional step actually converge faster than bisection?
|
|
|
|
The principle under test (Valeri's formulation):
|
|
in bisection the next point is the MIDPOINT of the bracket;
|
|
in the similarity version the next point is placed proportionally to how
|
|
SIMILAR the current result f(x) is to the expectation (zero):
|
|
similar -> small step
|
|
dissimilar-> large step
|
|
|
|
Formalisation used here. Bracket [a,b] with f(a)*f(b) < 0. Define the
|
|
similarity of an endpoint's value to the target 0 in a SCALE-FREE way, using
|
|
only the data inside the bracket (no parameter):
|
|
|
|
S(a) = |f(b)| / (|f(a)| + |f(b)|) (a is similar to the root when |f(a)| is small)
|
|
S(b) = |f(a)| / (|f(a)| + |f(b)|) S(a) + S(b) = 1
|
|
|
|
Next point = similarity-weighted aggregation of the prototypes {a, b}:
|
|
|
|
x = S(a)*a + S(b)*b
|
|
|
|
which is identical to stepping away from `a` by (1 - S(a)) of the bracket.
|
|
Bisection is the degenerate case S(a) = S(b) = 1/2, i.e. "no similarity
|
|
information used at all".
|
|
|
|
This file measures function evaluations to reach |x - root| <= tol, for
|
|
- bisection (S == 1/2, the null model)
|
|
- sim-step (the principle, raw)
|
|
and prints the per-problem numbers. No aggregate statistic is used to decide
|
|
anything; every problem is shown.
|
|
"""
|
|
from __future__ import annotations
|
|
|
|
import math
|
|
from dataclasses import dataclass
|
|
from typing import Callable
|
|
|
|
TOL = 1e-12
|
|
MAXFEV = 200
|
|
|
|
|
|
@dataclass
|
|
class Problem:
|
|
name: str
|
|
f: Callable[[float], float]
|
|
a: float
|
|
b: float
|
|
root: float
|
|
|
|
|
|
PROBLEMS = [
|
|
Problem("x^2 - 2 ", lambda x: x * x - 2.0, 0.0, 2.0, math.sqrt(2.0)),
|
|
Problem("cos(x) - x ", lambda x: math.cos(x) - x, 0.0, 1.0, 0.7390851332151607),
|
|
Problem("exp(x) - 2 ", lambda x: math.exp(x) - 2.0, 0.0, 2.0, math.log(2.0)),
|
|
Problem("x^3 - x - 2 ", lambda x: x ** 3 - x - 2.0, 1.0, 2.0, 1.5213797068045676),
|
|
Problem("sin(x) - 0.5 ", lambda x: math.sin(x) - 0.5, 0.0, 1.5, math.asin(0.5)),
|
|
# the classical stalling case for regula falsi: strong convexity, root near
|
|
# one end -> one endpoint never moves
|
|
Problem("x^10 - 1 (STALL case) ", lambda x: x ** 10 - 1.0, 0.0, 1.3, 1.0),
|
|
Problem("exp(-x)*(x-1) (flat end)", lambda x: math.exp(-x) * (x - 1.0), -1.0, 40.0, 1.0),
|
|
Problem("atan(x) ", lambda x: math.atan(x), -1.0, 3.0, 0.0),
|
|
]
|
|
|
|
|
|
def _count(f):
|
|
"""Wrap f with an evaluation counter."""
|
|
state = {"n": 0}
|
|
|
|
def g(x):
|
|
state["n"] += 1
|
|
return f(x)
|
|
|
|
return g, state
|
|
|
|
|
|
def bisection(f, a, b, tol=TOL, maxfev=MAXFEV):
|
|
g, st = _count(f)
|
|
fa, fb = g(a), g(b)
|
|
if fa == 0.0:
|
|
return a, st["n"]
|
|
if fb == 0.0:
|
|
return b, st["n"]
|
|
x = 0.5 * (a + b)
|
|
while st["n"] < maxfev:
|
|
x = 0.5 * (a + b)
|
|
fx = g(x)
|
|
if fx == 0.0 or (b - a) <= tol:
|
|
return x, st["n"]
|
|
if fa * fx < 0.0:
|
|
b, fb = x, fx
|
|
else:
|
|
a, fa = x, fx
|
|
return x, st["n"]
|
|
|
|
|
|
def sim_step(f, a, b, tol=TOL, maxfev=MAXFEV):
|
|
"""The raw principle: x = S(a)*a + S(b)*b with S scale-free inside the bracket."""
|
|
g, st = _count(f)
|
|
fa, fb = g(a), g(b)
|
|
if fa == 0.0:
|
|
return a, st["n"]
|
|
if fb == 0.0:
|
|
return b, st["n"]
|
|
x = 0.5 * (a + b)
|
|
while st["n"] < maxfev:
|
|
da, db = abs(fa), abs(fb)
|
|
tot = da + db
|
|
if tot == 0.0:
|
|
return 0.5 * (a + b), st["n"]
|
|
s_a = db / tot # similarity of endpoint a to the root
|
|
s_b = da / tot
|
|
x_new = s_a * a + s_b * b
|
|
# guard: keep the iterate strictly inside the bracket (numerical only,
|
|
# not a tuning parameter)
|
|
lo, hi = (a, b) if a < b else (b, a)
|
|
if not (lo < x_new < hi):
|
|
x_new = 0.5 * (a + b)
|
|
x = x_new
|
|
fx = g(x)
|
|
if fx == 0.0 or abs(b - a) <= tol:
|
|
return x, st["n"]
|
|
if fa * fx < 0.0:
|
|
b, fb = x, fx
|
|
else:
|
|
a, fa = x, fx
|
|
return x, st["n"]
|
|
|
|
|
|
def run():
|
|
print(f"{'problem':26s} {'bisect fev':>11s} {'bisect err':>11s} "
|
|
f"{'sim fev':>9s} {'sim err':>11s} verdict")
|
|
print("-" * 90)
|
|
wins = ties = losses = 0
|
|
for p in PROBLEMS:
|
|
xb, nb = bisection(p.f, p.a, p.b)
|
|
xs, ns = sim_step(p.f, p.a, p.b)
|
|
eb = abs(xb - p.root)
|
|
es = abs(xs - p.root)
|
|
# a method only "wins" if it reached the root at all
|
|
ok_b = eb <= 1e-8
|
|
ok_s = es <= 1e-8
|
|
if ok_s and (not ok_b or ns < nb):
|
|
verdict = "SIM faster"
|
|
wins += 1
|
|
elif ok_s and ok_b and ns == nb:
|
|
verdict = "tie"
|
|
ties += 1
|
|
else:
|
|
verdict = "SIM slower/failed"
|
|
losses += 1
|
|
print(f"{p.name:26s} {nb:11d} {eb:11.2e} {ns:9d} {es:11.2e} {verdict}")
|
|
print("-" * 90)
|
|
print(f"SIM faster: {wins} tie: {ties} SIM slower/failed: {losses}"
|
|
f" (of {len(PROBLEMS)}; maxfev={MAXFEV})")
|
|
|
|
|
|
if __name__ == "__main__":
|
|
run()
|