test_rootfinding.py 3.1 KB

12345678910111213141516171819202122232425262728293031323334353637383940414243444546474849505152535455565758596061626364656667686970717273747576777879808182838485868788899091
  1. import pytest
  2. from mpmath import *
  3. from mpmath.calculus.optimization import Secant, Muller, Bisection, Illinois, \
  4. Pegasus, Anderson, Ridder, ANewton, Newton, MNewton, MDNewton
  5. def test_findroot():
  6. # old tests, assuming secant
  7. mp.dps = 15
  8. assert findroot(lambda x: 4*x-3, mpf(5)).ae(0.75)
  9. assert findroot(sin, mpf(3)).ae(pi)
  10. assert findroot(sin, (mpf(3), mpf(3.14))).ae(pi)
  11. assert findroot(lambda x: x*x+1, mpc(2+2j)).ae(1j)
  12. # test all solvers with 1 starting point
  13. f = lambda x: cos(x)
  14. for solver in [Newton, Secant, MNewton, Muller, ANewton]:
  15. x = findroot(f, 2., solver=solver)
  16. assert abs(f(x)) < eps
  17. # test all solvers with interval of 2 points
  18. for solver in [Secant, Muller, Bisection, Illinois, Pegasus, Anderson,
  19. Ridder]:
  20. x = findroot(f, (1., 2.), solver=solver)
  21. assert abs(f(x)) < eps
  22. # test types
  23. f = lambda x: (x - 2)**2
  24. assert isinstance(findroot(f, 1, tol=1e-10), mpf)
  25. assert isinstance(iv.findroot(f, 1., tol=1e-10), iv.mpf)
  26. assert isinstance(fp.findroot(f, 1, tol=1e-10), float)
  27. assert isinstance(fp.findroot(f, 1+0j, tol=1e-10), complex)
  28. # issue 401
  29. with pytest.raises(ValueError):
  30. with workprec(2):
  31. findroot(lambda x: x**2 - 4456178*x + 60372201703370,
  32. mpc(real='5.278e+13', imag='-5.278e+13'))
  33. # issue 192
  34. with pytest.raises(ValueError):
  35. findroot(lambda x: -1, 0)
  36. # issue 387
  37. with pytest.raises(ValueError):
  38. findroot(lambda p: (1 - p)**30 - 1, 0.9)
  39. def test_bisection():
  40. # issue 273
  41. assert findroot(lambda x: x**2-1,(0,2),solver='bisect') == 1
  42. def test_mnewton():
  43. f = lambda x: polyval([1,3,3,1],x)
  44. x = findroot(f, -0.9, solver='mnewton')
  45. assert abs(f(x)) < eps
  46. def test_anewton():
  47. f = lambda x: (x - 2)**100
  48. x = findroot(f, 1., solver=ANewton)
  49. assert abs(f(x)) < eps
  50. def test_muller():
  51. f = lambda x: (2 + x)**3 + 2
  52. x = findroot(f, 1., solver=Muller)
  53. assert abs(f(x)) < eps
  54. def test_multiplicity():
  55. for i in range(1, 5):
  56. assert multiplicity(lambda x: (x - 1)**i, 1) == i
  57. assert multiplicity(lambda x: x**2, 1) == 0
  58. def test_multidimensional():
  59. def f(*x):
  60. return [3*x[0]**2-2*x[1]**2-1, x[0]**2-2*x[0]+x[1]**2+2*x[1]-8]
  61. assert mnorm(jacobian(f, (1,-2)) - matrix([[6,8],[0,-2]]),1) < 1.e-7
  62. for x, error in MDNewton(mp, f, (1,-2), verbose=0,
  63. norm=lambda x: norm(x, inf)):
  64. pass
  65. assert norm(f(*x), 2) < 1e-14
  66. # The Chinese mathematician Zhu Shijie was the very first to solve this
  67. # nonlinear system 700 years ago
  68. f1 = lambda x, y: -x + 2*y
  69. f2 = lambda x, y: (x**2 + x*(y**2 - 2) - 4*y) / (x + 4)
  70. f3 = lambda x, y: sqrt(x**2 + y**2)
  71. def f(x, y):
  72. f1x = f1(x, y)
  73. return (f2(x, y) - f1x, f3(x, y) - f1x)
  74. x = findroot(f, (10, 10))
  75. assert [int(round(i)) for i in x] == [3, 4]
  76. def test_trivial():
  77. assert findroot(lambda x: 0, 1) == 1
  78. assert findroot(lambda x: x, 0) == 0
  79. #assert findroot(lambda x, y: x + y, (1, -1)) == (1, -1)