test_elliptic_integrals.py 6.7 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179
  1. from sympy.core.numbers import (I, Rational, oo, pi, zoo)
  2. from sympy.core.singleton import S
  3. from sympy.core.symbol import (Dummy, Symbol)
  4. from sympy.functions.elementary.hyperbolic import atanh
  5. from sympy.functions.elementary.miscellaneous import sqrt
  6. from sympy.functions.elementary.trigonometric import (sin, tan)
  7. from sympy.functions.special.gamma_functions import gamma
  8. from sympy.functions.special.hyper import (hyper, meijerg)
  9. from sympy.integrals.integrals import Integral
  10. from sympy.series.order import O
  11. from sympy.functions.special.elliptic_integrals import (elliptic_k as K,
  12. elliptic_f as F, elliptic_e as E, elliptic_pi as P)
  13. from sympy.core.random import (test_derivative_numerically as td,
  14. random_complex_number as randcplx,
  15. verify_numerically as tn)
  16. from sympy.abc import z, m, n
  17. i = Symbol('i', integer=True)
  18. j = Symbol('k', integer=True, positive=True)
  19. t = Dummy('t')
  20. def test_K():
  21. assert K(0) == pi/2
  22. assert K(S.Half) == 8*pi**Rational(3, 2)/gamma(Rational(-1, 4))**2
  23. assert K(1) is zoo
  24. assert K(-1) == gamma(Rational(1, 4))**2/(4*sqrt(2*pi))
  25. assert K(oo) == 0
  26. assert K(-oo) == 0
  27. assert K(I*oo) == 0
  28. assert K(-I*oo) == 0
  29. assert K(zoo) == 0
  30. assert K(z).diff(z) == (E(z) - (1 - z)*K(z))/(2*z*(1 - z))
  31. assert td(K(z), z)
  32. zi = Symbol('z', real=False)
  33. assert K(zi).conjugate() == K(zi.conjugate())
  34. zr = Symbol('z', negative=True)
  35. assert K(zr).conjugate() == K(zr)
  36. assert K(z).rewrite(hyper) == \
  37. (pi/2)*hyper((S.Half, S.Half), (S.One,), z)
  38. assert tn(K(z), (pi/2)*hyper((S.Half, S.Half), (S.One,), z))
  39. assert K(z).rewrite(meijerg) == \
  40. meijerg(((S.Half, S.Half), []), ((S.Zero,), (S.Zero,)), -z)/2
  41. assert tn(K(z), meijerg(((S.Half, S.Half), []), ((S.Zero,), (S.Zero,)), -z)/2)
  42. assert K(z).series(z) == pi/2 + pi*z/8 + 9*pi*z**2/128 + \
  43. 25*pi*z**3/512 + 1225*pi*z**4/32768 + 3969*pi*z**5/131072 + O(z**6)
  44. assert K(m).rewrite(Integral).dummy_eq(
  45. Integral(1/sqrt(1 - m*sin(t)**2), (t, 0, pi/2)))
  46. def test_F():
  47. assert F(z, 0) == z
  48. assert F(0, m) == 0
  49. assert F(pi*i/2, m) == i*K(m)
  50. assert F(z, oo) == 0
  51. assert F(z, -oo) == 0
  52. assert F(-z, m) == -F(z, m)
  53. assert F(z, m).diff(z) == 1/sqrt(1 - m*sin(z)**2)
  54. assert F(z, m).diff(m) == E(z, m)/(2*m*(1 - m)) - F(z, m)/(2*m) - \
  55. sin(2*z)/(4*(1 - m)*sqrt(1 - m*sin(z)**2))
  56. r = randcplx()
  57. assert td(F(z, r), z)
  58. assert td(F(r, m), m)
  59. mi = Symbol('m', real=False)
  60. assert F(z, mi).conjugate() == F(z.conjugate(), mi.conjugate())
  61. mr = Symbol('m', negative=True)
  62. assert F(z, mr).conjugate() == F(z.conjugate(), mr)
  63. assert F(z, m).series(z) == \
  64. z + z**5*(3*m**2/40 - m/30) + m*z**3/6 + O(z**6)
  65. assert F(z, m).rewrite(Integral).dummy_eq(
  66. Integral(1/sqrt(1 - m*sin(t)**2), (t, 0, z)))
  67. def test_E():
  68. assert E(z, 0) == z
  69. assert E(0, m) == 0
  70. assert E(i*pi/2, m) == i*E(m)
  71. assert E(z, oo) is zoo
  72. assert E(z, -oo) is zoo
  73. assert E(0) == pi/2
  74. assert E(1) == 1
  75. assert E(oo) == I*oo
  76. assert E(-oo) is oo
  77. assert E(zoo) is zoo
  78. assert E(-z, m) == -E(z, m)
  79. assert E(z, m).diff(z) == sqrt(1 - m*sin(z)**2)
  80. assert E(z, m).diff(m) == (E(z, m) - F(z, m))/(2*m)
  81. assert E(z).diff(z) == (E(z) - K(z))/(2*z)
  82. r = randcplx()
  83. assert td(E(r, m), m)
  84. assert td(E(z, r), z)
  85. assert td(E(z), z)
  86. mi = Symbol('m', real=False)
  87. assert E(z, mi).conjugate() == E(z.conjugate(), mi.conjugate())
  88. assert E(mi).conjugate() == E(mi.conjugate())
  89. mr = Symbol('m', negative=True)
  90. assert E(z, mr).conjugate() == E(z.conjugate(), mr)
  91. assert E(mr).conjugate() == E(mr)
  92. assert E(z).rewrite(hyper) == (pi/2)*hyper((Rational(-1, 2), S.Half), (S.One,), z)
  93. assert tn(E(z), (pi/2)*hyper((Rational(-1, 2), S.Half), (S.One,), z))
  94. assert E(z).rewrite(meijerg) == \
  95. -meijerg(((S.Half, Rational(3, 2)), []), ((S.Zero,), (S.Zero,)), -z)/4
  96. assert tn(E(z), -meijerg(((S.Half, Rational(3, 2)), []), ((S.Zero,), (S.Zero,)), -z)/4)
  97. assert E(z, m).series(z) == \
  98. z + z**5*(-m**2/40 + m/30) - m*z**3/6 + O(z**6)
  99. assert E(z).series(z) == pi/2 - pi*z/8 - 3*pi*z**2/128 - \
  100. 5*pi*z**3/512 - 175*pi*z**4/32768 - 441*pi*z**5/131072 + O(z**6)
  101. assert E(z, m).rewrite(Integral).dummy_eq(
  102. Integral(sqrt(1 - m*sin(t)**2), (t, 0, z)))
  103. assert E(m).rewrite(Integral).dummy_eq(
  104. Integral(sqrt(1 - m*sin(t)**2), (t, 0, pi/2)))
  105. def test_P():
  106. assert P(0, z, m) == F(z, m)
  107. assert P(1, z, m) == F(z, m) + \
  108. (sqrt(1 - m*sin(z)**2)*tan(z) - E(z, m))/(1 - m)
  109. assert P(n, i*pi/2, m) == i*P(n, m)
  110. assert P(n, z, 0) == atanh(sqrt(n - 1)*tan(z))/sqrt(n - 1)
  111. assert P(n, z, n) == F(z, n) - P(1, z, n) + tan(z)/sqrt(1 - n*sin(z)**2)
  112. assert P(oo, z, m) == 0
  113. assert P(-oo, z, m) == 0
  114. assert P(n, z, oo) == 0
  115. assert P(n, z, -oo) == 0
  116. assert P(0, m) == K(m)
  117. assert P(1, m) is zoo
  118. assert P(n, 0) == pi/(2*sqrt(1 - n))
  119. assert P(2, 1) is -oo
  120. assert P(-1, 1) is oo
  121. assert P(n, n) == E(n)/(1 - n)
  122. assert P(n, -z, m) == -P(n, z, m)
  123. ni, mi = Symbol('n', real=False), Symbol('m', real=False)
  124. assert P(ni, z, mi).conjugate() == \
  125. P(ni.conjugate(), z.conjugate(), mi.conjugate())
  126. nr, mr = Symbol('n', negative=True), \
  127. Symbol('m', negative=True)
  128. assert P(nr, z, mr).conjugate() == P(nr, z.conjugate(), mr)
  129. assert P(n, m).conjugate() == P(n.conjugate(), m.conjugate())
  130. assert P(n, z, m).diff(n) == (E(z, m) + (m - n)*F(z, m)/n +
  131. (n**2 - m)*P(n, z, m)/n - n*sqrt(1 -
  132. m*sin(z)**2)*sin(2*z)/(2*(1 - n*sin(z)**2)))/(2*(m - n)*(n - 1))
  133. assert P(n, z, m).diff(z) == 1/(sqrt(1 - m*sin(z)**2)*(1 - n*sin(z)**2))
  134. assert P(n, z, m).diff(m) == (E(z, m)/(m - 1) + P(n, z, m) -
  135. m*sin(2*z)/(2*(m - 1)*sqrt(1 - m*sin(z)**2)))/(2*(n - m))
  136. assert P(n, m).diff(n) == (E(m) + (m - n)*K(m)/n +
  137. (n**2 - m)*P(n, m)/n)/(2*(m - n)*(n - 1))
  138. assert P(n, m).diff(m) == (E(m)/(m - 1) + P(n, m))/(2*(n - m))
  139. # These tests fail due to
  140. # https://github.com/fredrik-johansson/mpmath/issues/571#issuecomment-777201962
  141. # https://github.com/sympy/sympy/issues/20933#issuecomment-777080385
  142. #
  143. # rx, ry = randcplx(), randcplx()
  144. # assert td(P(n, rx, ry), n)
  145. # assert td(P(rx, z, ry), z)
  146. # assert td(P(rx, ry, m), m)
  147. assert P(n, z, m).series(z) == z + z**3*(m/6 + n/3) + \
  148. z**5*(3*m**2/40 + m*n/10 - m/30 + n**2/5 - n/15) + O(z**6)
  149. assert P(n, z, m).rewrite(Integral).dummy_eq(
  150. Integral(1/((1 - n*sin(t)**2)*sqrt(1 - m*sin(t)**2)), (t, 0, z)))
  151. assert P(n, m).rewrite(Integral).dummy_eq(
  152. Integral(1/((1 - n*sin(t)**2)*sqrt(1 - m*sin(t)**2)), (t, 0, pi/2)))