Files
2026-07-13 12:32:53 +08:00

177 lines
6.2 KiB
Python

import pytest
from mpmath import (cos, eps, findroot, fp, inf, iv, jacobian, matrix, mnorm,
mp, mpc, mpf, multiplicity, norm, pi, polyval, sin, sqrt,
workprec)
from mpmath.calculus.optimization import (Anderson, ANewton, Bisection,
Illinois, MDNewton, MNewton, Muller,
Newton, Pegasus, Ridder, Secant, ModAB, Brent)
def test_findroot():
# old tests, assuming secant
assert findroot(lambda x: 4*x-3, mpf(5)).ae(0.75)
assert findroot(sin, mpf(3)).ae(pi)
assert findroot(sin, (mpf(3), mpf(3.14))).ae(pi)
assert findroot(lambda x: x*x+1, mpc(2+2j)).ae(1j)
# test all solvers with 1 starting point
f = lambda x: cos(x)
for solver in [Newton, Secant, MNewton, Muller, ANewton]:
x = findroot(f, 2., solver=solver)
assert abs(f(x)) < eps
# test all solvers with interval of 2 points
for solver in [Secant, Muller, Bisection, Illinois, Pegasus, Anderson,
Ridder, ModAB, Brent]:
x = findroot(f, (1., 2.), solver=solver)
assert abs(f(x)) < eps
# test types
f = lambda x: (x - 2)**2
assert isinstance(findroot(f, 1, tol=1e-10), mpf)
assert isinstance(iv.findroot(f, 1., tol=1e-10), iv.mpf)
assert isinstance(fp.findroot(f, 1, tol=1e-10), float)
assert isinstance(fp.findroot(f, 1+0j, tol=1e-10), complex)
# issue 401
with pytest.raises(ValueError):
with workprec(2):
findroot(lambda x: x**2 - 4456178*x + 60372201703370,
mpc(real='5.278e+13', imag='-5.278e+13'))
# issue 192
with pytest.raises(ValueError):
findroot(lambda x: -1, 0)
# issue 387
with pytest.raises(ValueError):
findroot(lambda p: (1 - p)**30 - 1, 0.9)
def test_bisection():
# issue 273
assert findroot(lambda x: x**2-1,(0,2),solver='bisect') == 1
with pytest.raises(ValueError):
findroot(lambda x: x**2-1, (4, 2), solver='bisect') == 1
# issue 285
mp.dps = 240
sol = -mp.ceil(mp.log(abs(findroot(lambda x: mp.sign(x - 3), (1, 4),
solver='bisect', verify=False,
tol=1e-200) - 3))/mp.log(10))
assert sol.ae(200)
# issue 339
mp.dps = 15
res = mpf('0.73908513321516064')
for dps in [100, 200, 300, 1000]:
with mp.workdps(dps):
sol = findroot(lambda x: cos(x) - x, [0, 1], solver='bisect')
assert (+sol).ae(res)
def test_mnewton():
f = lambda x: polyval([1, 3, 3, 1], x)
x = findroot(f, -0.9, solver='mnewton')
assert abs(f(x)) < eps
def test_anewton():
f = lambda x: (x - 2)**100
x = findroot(f, 1., solver=ANewton)
assert abs(f(x)) < eps
def test_muller():
f = lambda x: (2 + x)**3 + 2
x = findroot(f, 1., solver=Muller)
assert abs(f(x)) < eps
def test_ridder():
f = lambda x: cos(x)/x
x = findroot(f, (1, 2), solver='ridder')
assert abs(f(x)) < eps
def test_brent():
f = lambda x: cos(x)/x
x = findroot(f, (1, 2), solver='brent')
assert abs(f(x)) < eps
with pytest.raises(ValueError, match="expected interval of 2 points"):
findroot(lambda x: x**2 - 1, (0,), solver='brent')
with pytest.raises(ValueError, match="Function must have opposite signs"):
findroot(lambda x: x**2 - 1, (2, 4), solver='brent')
assert findroot(lambda x: x, (-1, 2), solver='brent') == 0.0
assert findroot(lambda x: x, (-1, 1), solver='brent') == 0.0
def test_modAB():
assert findroot(lambda x: x**2 - 1, (0, 2), solver='modAB') == 1
# test ordering
assert findroot(lambda x: x**2 - 1, (2, 0), solver='modAB') == 1
with pytest.raises(ValueError, match="expected interval of 2 points"):
findroot(lambda x: x**2 - 1, (0,), solver='modAB')
with pytest.raises(ValueError, match="Function must have opposite signs"):
findroot(lambda x: x**2 - 1, (2, 4), solver='modAB')
# test exact zero hit
assert findroot(lambda x: x, (-1, 1), solver='modAB') == 0.0
# test bisection to secant switch for a purely linear function
f_linear = lambda x: 2*x - 4
assert mp.almosteq(findroot(f_linear, (0, 5), solver='modAB'), 2.0)
f_convex = lambda x: x**10 - 1
assert mp.almosteq(findroot(f_convex, (0.1, 2.0), solver='modAB'), 1.0)
f_concave = lambda x: 1 - x**10
assert mp.almosteq(findroot(f_concave, (2.0, 0.1), solver='modAB'), 1.0)
f_cubic_inflection = lambda x: x**3 - 3*x + 3
root = findroot(f_cubic_inflection, (-3, 2), solver='modAB')
assert abs(f_cubic_inflection(root)) < eps
# test reset to Bisection if the interval width exceeds the threshold
f_step = lambda x: mp.sin(x) if x > 1 else x - 1
assert mp.almosteq(findroot(f_step, (0.4, 3.0), solver='modAB'), 1.0)
def test_multiplicity():
for i in range(1, 5):
assert multiplicity(lambda x: (x - 1)**i, 1) == i
assert multiplicity(lambda x: x**2, 1) == 0
def test_multidimensional(capsys):
def f(*x):
return [3*x[0]**2-2*x[1]**2-1, x[0]**2-2*x[0]+x[1]**2+2*x[1]-8]
assert mnorm(jacobian(f, (1,-2)) - matrix([[6,8],[0,-2]]),1) < 1.e-7
for x, error in MDNewton(mp, f, (1,-2), verbose=0,
norm=lambda x: norm(x, inf)):
pass
assert norm(f(*x), 2) < 1e-14
for x, error in MDNewton(mp, f, (1,-2), verbose=1,
norm=lambda x: norm(x, inf)):
pass
assert norm(f(*x), 2) < 1e-14
captured = capsys.readouterr()
assert captured.out.find("canceled, won't get more exact") >= 0
# The Chinese mathematician Zhu Shijie was the very first to solve this
# nonlinear system 700 years ago
f1 = lambda x, y: -x + 2*y
f2 = lambda x, y: (x**2 + x*(y**2 - 2) - 4*y) / (x + 4)
f3 = lambda x, y: sqrt(x**2 + y**2)
def f(x, y):
f1x = f1(x, y)
return (f2(x, y) - f1x, f3(x, y) - f1x)
x = findroot(f, (10, 10))
assert [round(i) for i in x] == [3, 4]
def test_trivial():
assert findroot(lambda x: 0, 1) == 1
assert findroot(lambda x: x, 0) == 0
#assert findroot(lambda x, y: x + y, (1, -1)) == (1, -1)
def test_issue_869():
f = [lambda x: sqrt(x) + 1]
pytest.raises(mp.ComplexResult, lambda: findroot(f, [-1]))