test_quad.py 3.8 KB

1234567891011121314151617181920212223242526272829303132333435363738394041424344454647484950515253545556575859606162636465666768697071727374757677787980818283848586878889909192939495
  1. import pytest
  2. from mpmath import *
  3. def ae(a, b):
  4. return abs(a-b) < 10**(-mp.dps+5)
  5. def test_basic_integrals():
  6. for prec in [15, 30, 100]:
  7. mp.dps = prec
  8. assert ae(quadts(lambda x: x**3 - 3*x**2, [-2, 4]), -12)
  9. assert ae(quadgl(lambda x: x**3 - 3*x**2, [-2, 4]), -12)
  10. assert ae(quadts(sin, [0, pi]), 2)
  11. assert ae(quadts(sin, [0, 2*pi]), 0)
  12. assert ae(quadts(exp, [-inf, -1]), 1/e)
  13. assert ae(quadts(lambda x: exp(-x), [0, inf]), 1)
  14. assert ae(quadts(lambda x: exp(-x*x), [-inf, inf]), sqrt(pi))
  15. assert ae(quadts(lambda x: 1/(1+x*x), [-1, 1]), pi/2)
  16. assert ae(quadts(lambda x: 1/(1+x*x), [-inf, inf]), pi)
  17. assert ae(quadts(lambda x: 2*sqrt(1-x*x), [-1, 1]), pi)
  18. mp.dps = 15
  19. def test_multiple_intervals():
  20. y,err = quad(lambda x: sign(x), [-0.5, 0.9, 1], maxdegree=2, error=True)
  21. assert abs(y-0.5) < 2*err
  22. def test_quad_symmetry():
  23. assert quadts(sin, [-1, 1]) == 0
  24. assert quadgl(sin, [-1, 1]) == 0
  25. def test_quad_infinite_mirror():
  26. # Check mirrored infinite interval
  27. assert ae(quad(lambda x: exp(-x*x), [inf,-inf]), -sqrt(pi))
  28. assert ae(quad(lambda x: exp(x), [0,-inf]), -1)
  29. def test_quadgl_linear():
  30. assert quadgl(lambda x: x, [0, 1], maxdegree=1).ae(0.5)
  31. def test_complex_integration():
  32. assert quadts(lambda x: x, [0, 1+j]).ae(j)
  33. def test_quadosc():
  34. mp.dps = 15
  35. assert quadosc(lambda x: sin(x)/x, [0, inf], period=2*pi).ae(pi/2)
  36. # Double integrals
  37. def test_double_trivial():
  38. assert ae(quadts(lambda x, y: x, [0, 1], [0, 1]), 0.5)
  39. assert ae(quadts(lambda x, y: x, [-1, 1], [-1, 1]), 0.0)
  40. def test_double_1():
  41. assert ae(quadts(lambda x, y: cos(x+y/2), [-pi/2, pi/2], [0, pi]), 4)
  42. def test_double_2():
  43. assert ae(quadts(lambda x, y: (x-1)/((1-x*y)*log(x*y)), [0, 1], [0, 1]), euler)
  44. def test_double_3():
  45. assert ae(quadts(lambda x, y: 1/sqrt(1+x*x+y*y), [-1, 1], [-1, 1]), 4*log(2+sqrt(3))-2*pi/3)
  46. def test_double_4():
  47. assert ae(quadts(lambda x, y: 1/(1-x*x * y*y), [0, 1], [0, 1]), pi**2 / 8)
  48. def test_double_5():
  49. assert ae(quadts(lambda x, y: 1/(1-x*y), [0, 1], [0, 1]), pi**2 / 6)
  50. def test_double_6():
  51. assert ae(quadts(lambda x, y: exp(-(x+y)), [0, inf], [0, inf]), 1)
  52. def test_double_7():
  53. assert ae(quadts(lambda x, y: exp(-x*x-y*y), [-inf, inf], [-inf, inf]), pi)
  54. # Test integrals from "Experimentation in Mathematics" by Borwein,
  55. # Bailey & Girgensohn
  56. def test_expmath_integrals():
  57. for prec in [15, 30, 50]:
  58. mp.dps = prec
  59. assert ae(quadts(lambda x: x/sinh(x), [0, inf]), pi**2 / 4)
  60. assert ae(quadts(lambda x: log(x)**2 / (1+x**2), [0, inf]), pi**3 / 8)
  61. assert ae(quadts(lambda x: (1+x**2)/(1+x**4), [0, inf]), pi/sqrt(2))
  62. assert ae(quadts(lambda x: log(x)/cosh(x)**2, [0, inf]), log(pi)-2*log(2)-euler)
  63. assert ae(quadts(lambda x: log(1+x**3)/(1-x+x**2), [0, inf]), 2*pi*log(3)/sqrt(3))
  64. assert ae(quadts(lambda x: log(x)**2 / (x**2+x+1), [0, 1]), 8*pi**3 / (81*sqrt(3)))
  65. assert ae(quadts(lambda x: log(cos(x))**2, [0, pi/2]), pi/2 * (log(2)**2+pi**2/12))
  66. assert ae(quadts(lambda x: x**2 / sin(x)**2, [0, pi/2]), pi*log(2))
  67. assert ae(quadts(lambda x: x**2/sqrt(exp(x)-1), [0, inf]), 4*pi*(log(2)**2 + pi**2/12))
  68. assert ae(quadts(lambda x: x*exp(-x)*sqrt(1-exp(-2*x)), [0, inf]), pi*(1+2*log(2))/8)
  69. mp.dps = 15
  70. # Do not reach full accuracy
  71. @pytest.mark.xfail
  72. def test_expmath_fail():
  73. assert ae(quadts(lambda x: sqrt(tan(x)), [0, pi/2]), pi*sqrt(2)/2)
  74. assert ae(quadts(lambda x: atan(x)/(x*sqrt(1-x**2)), [0, 1]), pi*log(1+sqrt(2))/2)
  75. assert ae(quadts(lambda x: log(1+x**2)/x**2, [0, 1]), pi/2-log(2))
  76. assert ae(quadts(lambda x: x**2/((1+x**4)*sqrt(1-x**4)), [0, 1]), pi/8)