elliptic_integrals.py 15 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445
  1. """ Elliptic Integrals. """
  2. from sympy.core import S, pi, I, Rational
  3. from sympy.core.function import DefinedFunction, ArgumentIndexError
  4. from sympy.core.symbol import Dummy,uniquely_named_symbol
  5. from sympy.functions.elementary.complexes import sign
  6. from sympy.functions.elementary.hyperbolic import atanh
  7. from sympy.functions.elementary.miscellaneous import sqrt
  8. from sympy.functions.elementary.trigonometric import sin, tan
  9. from sympy.functions.special.gamma_functions import gamma
  10. from sympy.functions.special.hyper import hyper, meijerg
  11. class elliptic_k(DefinedFunction):
  12. r"""
  13. The complete elliptic integral of the first kind, defined by
  14. .. math:: K(m) = F\left(\tfrac{\pi}{2}\middle| m\right)
  15. where $F\left(z\middle| m\right)$ is the Legendre incomplete
  16. elliptic integral of the first kind.
  17. Explanation
  18. ===========
  19. The function $K(m)$ is a single-valued function on the complex
  20. plane with branch cut along the interval $(1, \infty)$.
  21. Note that our notation defines the incomplete elliptic integral
  22. in terms of the parameter $m$ instead of the elliptic modulus
  23. (eccentricity) $k$.
  24. In this case, the parameter $m$ is defined as $m=k^2$.
  25. Examples
  26. ========
  27. >>> from sympy import elliptic_k, I
  28. >>> from sympy.abc import m
  29. >>> elliptic_k(0)
  30. pi/2
  31. >>> elliptic_k(1.0 + I)
  32. 1.50923695405127 + 0.625146415202697*I
  33. >>> elliptic_k(m).series(n=3)
  34. pi/2 + pi*m/8 + 9*pi*m**2/128 + O(m**3)
  35. See Also
  36. ========
  37. elliptic_f
  38. References
  39. ==========
  40. .. [1] https://en.wikipedia.org/wiki/Elliptic_integrals
  41. .. [2] https://functions.wolfram.com/EllipticIntegrals/EllipticK
  42. """
  43. @classmethod
  44. def eval(cls, m):
  45. if m.is_zero:
  46. return pi*S.Half
  47. elif m is S.Half:
  48. return 8*pi**Rational(3, 2)/gamma(Rational(-1, 4))**2
  49. elif m is S.One:
  50. return S.ComplexInfinity
  51. elif m is S.NegativeOne:
  52. return gamma(Rational(1, 4))**2/(4*sqrt(2*pi))
  53. elif m in (S.Infinity, S.NegativeInfinity, I*S.Infinity,
  54. I*S.NegativeInfinity, S.ComplexInfinity):
  55. return S.Zero
  56. def fdiff(self, argindex=1):
  57. m = self.args[0]
  58. return (elliptic_e(m) - (1 - m)*elliptic_k(m))/(2*m*(1 - m))
  59. def _eval_conjugate(self):
  60. m = self.args[0]
  61. if (m.is_real and (m - 1).is_positive) is False:
  62. return self.func(m.conjugate())
  63. def _eval_nseries(self, x, n, logx, cdir=0):
  64. from sympy.simplify import hyperexpand
  65. return hyperexpand(self.rewrite(hyper)._eval_nseries(x, n=n, logx=logx))
  66. def _eval_rewrite_as_hyper(self, m, **kwargs):
  67. return pi*S.Half*hyper((S.Half, S.Half), (S.One,), m)
  68. def _eval_rewrite_as_meijerg(self, m, **kwargs):
  69. return meijerg(((S.Half, S.Half), []), ((S.Zero,), (S.Zero,)), -m)/2
  70. def _eval_is_zero(self):
  71. m = self.args[0]
  72. if m.is_infinite:
  73. return True
  74. def _eval_rewrite_as_Integral(self, *args, **kwargs):
  75. from sympy.integrals.integrals import Integral
  76. t = Dummy(uniquely_named_symbol('t', args).name)
  77. m = self.args[0]
  78. return Integral(1/sqrt(1 - m*sin(t)**2), (t, 0, pi/2))
  79. class elliptic_f(DefinedFunction):
  80. r"""
  81. The Legendre incomplete elliptic integral of the first
  82. kind, defined by
  83. .. math:: F\left(z\middle| m\right) =
  84. \int_0^z \frac{dt}{\sqrt{1 - m \sin^2 t}}
  85. Explanation
  86. ===========
  87. This function reduces to a complete elliptic integral of
  88. the first kind, $K(m)$, when $z = \pi/2$.
  89. Note that our notation defines the incomplete elliptic integral
  90. in terms of the parameter $m$ instead of the elliptic modulus
  91. (eccentricity) $k$.
  92. In this case, the parameter $m$ is defined as $m=k^2$.
  93. Examples
  94. ========
  95. >>> from sympy import elliptic_f, I
  96. >>> from sympy.abc import z, m
  97. >>> elliptic_f(z, m).series(z)
  98. z + z**5*(3*m**2/40 - m/30) + m*z**3/6 + O(z**6)
  99. >>> elliptic_f(3.0 + I/2, 1.0 + I)
  100. 2.909449841483 + 1.74720545502474*I
  101. See Also
  102. ========
  103. elliptic_k
  104. References
  105. ==========
  106. .. [1] https://en.wikipedia.org/wiki/Elliptic_integrals
  107. .. [2] https://functions.wolfram.com/EllipticIntegrals/EllipticF
  108. """
  109. @classmethod
  110. def eval(cls, z, m):
  111. if z.is_zero:
  112. return S.Zero
  113. if m.is_zero:
  114. return z
  115. k = 2*z/pi
  116. if k.is_integer:
  117. return k*elliptic_k(m)
  118. elif m in (S.Infinity, S.NegativeInfinity):
  119. return S.Zero
  120. elif z.could_extract_minus_sign():
  121. return -elliptic_f(-z, m)
  122. def fdiff(self, argindex=1):
  123. z, m = self.args
  124. fm = sqrt(1 - m*sin(z)**2)
  125. if argindex == 1:
  126. return 1/fm
  127. elif argindex == 2:
  128. return (elliptic_e(z, m)/(2*m*(1 - m)) - elliptic_f(z, m)/(2*m) -
  129. sin(2*z)/(4*(1 - m)*fm))
  130. raise ArgumentIndexError(self, argindex)
  131. def _eval_conjugate(self):
  132. z, m = self.args
  133. if (m.is_real and (m - 1).is_positive) is False:
  134. return self.func(z.conjugate(), m.conjugate())
  135. def _eval_rewrite_as_Integral(self, *args, **kwargs):
  136. from sympy.integrals.integrals import Integral
  137. t = Dummy(uniquely_named_symbol('t', args).name)
  138. z, m = self.args[0], self.args[1]
  139. return Integral(1/(sqrt(1 - m*sin(t)**2)), (t, 0, z))
  140. def _eval_is_zero(self):
  141. z, m = self.args
  142. if z.is_zero:
  143. return True
  144. if m.is_extended_real and m.is_infinite:
  145. return True
  146. class elliptic_e(DefinedFunction):
  147. r"""
  148. Called with two arguments $z$ and $m$, evaluates the
  149. incomplete elliptic integral of the second kind, defined by
  150. .. math:: E\left(z\middle| m\right) = \int_0^z \sqrt{1 - m \sin^2 t} dt
  151. Called with a single argument $m$, evaluates the Legendre complete
  152. elliptic integral of the second kind
  153. .. math:: E(m) = E\left(\tfrac{\pi}{2}\middle| m\right)
  154. Explanation
  155. ===========
  156. The function $E(m)$ is a single-valued function on the complex
  157. plane with branch cut along the interval $(1, \infty)$.
  158. Note that our notation defines the incomplete elliptic integral
  159. in terms of the parameter $m$ instead of the elliptic modulus
  160. (eccentricity) $k$.
  161. In this case, the parameter $m$ is defined as $m=k^2$.
  162. Examples
  163. ========
  164. >>> from sympy import elliptic_e, I
  165. >>> from sympy.abc import z, m
  166. >>> elliptic_e(z, m).series(z)
  167. z + z**5*(-m**2/40 + m/30) - m*z**3/6 + O(z**6)
  168. >>> elliptic_e(m).series(n=4)
  169. pi/2 - pi*m/8 - 3*pi*m**2/128 - 5*pi*m**3/512 + O(m**4)
  170. >>> elliptic_e(1 + I, 2 - I/2).n()
  171. 1.55203744279187 + 0.290764986058437*I
  172. >>> elliptic_e(0)
  173. pi/2
  174. >>> elliptic_e(2.0 - I)
  175. 0.991052601328069 + 0.81879421395609*I
  176. References
  177. ==========
  178. .. [1] https://en.wikipedia.org/wiki/Elliptic_integrals
  179. .. [2] https://functions.wolfram.com/EllipticIntegrals/EllipticE2
  180. .. [3] https://functions.wolfram.com/EllipticIntegrals/EllipticE
  181. """
  182. @classmethod
  183. def eval(cls, m, z=None):
  184. if z is not None:
  185. z, m = m, z
  186. k = 2*z/pi
  187. if m.is_zero:
  188. return z
  189. if z.is_zero:
  190. return S.Zero
  191. elif k.is_integer:
  192. return k*elliptic_e(m)
  193. elif m in (S.Infinity, S.NegativeInfinity):
  194. return S.ComplexInfinity
  195. elif z.could_extract_minus_sign():
  196. return -elliptic_e(-z, m)
  197. else:
  198. if m.is_zero:
  199. return pi/2
  200. elif m is S.One:
  201. return S.One
  202. elif m is S.Infinity:
  203. return I*S.Infinity
  204. elif m is S.NegativeInfinity:
  205. return S.Infinity
  206. elif m is S.ComplexInfinity:
  207. return S.ComplexInfinity
  208. def fdiff(self, argindex=1):
  209. if len(self.args) == 2:
  210. z, m = self.args
  211. if argindex == 1:
  212. return sqrt(1 - m*sin(z)**2)
  213. elif argindex == 2:
  214. return (elliptic_e(z, m) - elliptic_f(z, m))/(2*m)
  215. else:
  216. m = self.args[0]
  217. if argindex == 1:
  218. return (elliptic_e(m) - elliptic_k(m))/(2*m)
  219. raise ArgumentIndexError(self, argindex)
  220. def _eval_conjugate(self):
  221. if len(self.args) == 2:
  222. z, m = self.args
  223. if (m.is_real and (m - 1).is_positive) is False:
  224. return self.func(z.conjugate(), m.conjugate())
  225. else:
  226. m = self.args[0]
  227. if (m.is_real and (m - 1).is_positive) is False:
  228. return self.func(m.conjugate())
  229. def _eval_nseries(self, x, n, logx, cdir=0):
  230. from sympy.simplify import hyperexpand
  231. if len(self.args) == 1:
  232. return hyperexpand(self.rewrite(hyper)._eval_nseries(x, n=n, logx=logx))
  233. return super()._eval_nseries(x, n=n, logx=logx)
  234. def _eval_rewrite_as_hyper(self, *args, **kwargs):
  235. if len(args) == 1:
  236. m = args[0]
  237. return (pi/2)*hyper((Rational(-1, 2), S.Half), (S.One,), m)
  238. def _eval_rewrite_as_meijerg(self, *args, **kwargs):
  239. if len(args) == 1:
  240. m = args[0]
  241. return -meijerg(((S.Half, Rational(3, 2)), []), \
  242. ((S.Zero,), (S.Zero,)), -m)/4
  243. def _eval_rewrite_as_Integral(self, *args, **kwargs):
  244. from sympy.integrals.integrals import Integral
  245. z, m = (pi/2, self.args[0]) if len(self.args) == 1 else self.args
  246. t = Dummy(uniquely_named_symbol('t', args).name)
  247. return Integral(sqrt(1 - m*sin(t)**2), (t, 0, z))
  248. class elliptic_pi(DefinedFunction):
  249. r"""
  250. Called with three arguments $n$, $z$ and $m$, evaluates the
  251. Legendre incomplete elliptic integral of the third kind, defined by
  252. .. math:: \Pi\left(n; z\middle| m\right) = \int_0^z \frac{dt}
  253. {\left(1 - n \sin^2 t\right) \sqrt{1 - m \sin^2 t}}
  254. Called with two arguments $n$ and $m$, evaluates the complete
  255. elliptic integral of the third kind:
  256. .. math:: \Pi\left(n\middle| m\right) =
  257. \Pi\left(n; \tfrac{\pi}{2}\middle| m\right)
  258. Explanation
  259. ===========
  260. Note that our notation defines the incomplete elliptic integral
  261. in terms of the parameter $m$ instead of the elliptic modulus
  262. (eccentricity) $k$.
  263. In this case, the parameter $m$ is defined as $m=k^2$.
  264. Examples
  265. ========
  266. >>> from sympy import elliptic_pi, I
  267. >>> from sympy.abc import z, n, m
  268. >>> elliptic_pi(n, z, m).series(z, n=4)
  269. z + z**3*(m/6 + n/3) + O(z**4)
  270. >>> elliptic_pi(0.5 + I, 1.0 - I, 1.2)
  271. 2.50232379629182 - 0.760939574180767*I
  272. >>> elliptic_pi(0, 0)
  273. pi/2
  274. >>> elliptic_pi(1.0 - I/3, 2.0 + I)
  275. 3.29136443417283 + 0.32555634906645*I
  276. References
  277. ==========
  278. .. [1] https://en.wikipedia.org/wiki/Elliptic_integrals
  279. .. [2] https://functions.wolfram.com/EllipticIntegrals/EllipticPi3
  280. .. [3] https://functions.wolfram.com/EllipticIntegrals/EllipticPi
  281. """
  282. @classmethod
  283. def eval(cls, n, m, z=None):
  284. if z is not None:
  285. z, m = m, z
  286. if n.is_zero:
  287. return elliptic_f(z, m)
  288. elif n is S.One:
  289. return (elliptic_f(z, m) +
  290. (sqrt(1 - m*sin(z)**2)*tan(z) -
  291. elliptic_e(z, m))/(1 - m))
  292. k = 2*z/pi
  293. if k.is_integer:
  294. return k*elliptic_pi(n, m)
  295. elif m.is_zero:
  296. return atanh(sqrt(n - 1)*tan(z))/sqrt(n - 1)
  297. elif n == m:
  298. return (elliptic_f(z, n) - elliptic_pi(1, z, n) +
  299. tan(z)/sqrt(1 - n*sin(z)**2))
  300. elif n in (S.Infinity, S.NegativeInfinity):
  301. return S.Zero
  302. elif m in (S.Infinity, S.NegativeInfinity):
  303. return S.Zero
  304. elif z.could_extract_minus_sign():
  305. return -elliptic_pi(n, -z, m)
  306. if n.is_zero:
  307. return elliptic_f(z, m)
  308. if m.is_extended_real and m.is_infinite or \
  309. n.is_extended_real and n.is_infinite:
  310. return S.Zero
  311. else:
  312. if n.is_zero:
  313. return elliptic_k(m)
  314. elif n is S.One:
  315. return S.ComplexInfinity
  316. elif m.is_zero:
  317. return pi/(2*sqrt(1 - n))
  318. elif m == S.One:
  319. return S.NegativeInfinity/sign(n - 1)
  320. elif n == m:
  321. return elliptic_e(n)/(1 - n)
  322. elif n in (S.Infinity, S.NegativeInfinity):
  323. return S.Zero
  324. elif m in (S.Infinity, S.NegativeInfinity):
  325. return S.Zero
  326. if n.is_zero:
  327. return elliptic_k(m)
  328. if m.is_extended_real and m.is_infinite or \
  329. n.is_extended_real and n.is_infinite:
  330. return S.Zero
  331. def _eval_conjugate(self):
  332. if len(self.args) == 3:
  333. n, z, m = self.args
  334. if (n.is_real and (n - 1).is_positive) is False and \
  335. (m.is_real and (m - 1).is_positive) is False:
  336. return self.func(n.conjugate(), z.conjugate(), m.conjugate())
  337. else:
  338. n, m = self.args
  339. return self.func(n.conjugate(), m.conjugate())
  340. def fdiff(self, argindex=1):
  341. if len(self.args) == 3:
  342. n, z, m = self.args
  343. fm, fn = sqrt(1 - m*sin(z)**2), 1 - n*sin(z)**2
  344. if argindex == 1:
  345. return (elliptic_e(z, m) + (m - n)*elliptic_f(z, m)/n +
  346. (n**2 - m)*elliptic_pi(n, z, m)/n -
  347. n*fm*sin(2*z)/(2*fn))/(2*(m - n)*(n - 1))
  348. elif argindex == 2:
  349. return 1/(fm*fn)
  350. elif argindex == 3:
  351. return (elliptic_e(z, m)/(m - 1) +
  352. elliptic_pi(n, z, m) -
  353. m*sin(2*z)/(2*(m - 1)*fm))/(2*(n - m))
  354. else:
  355. n, m = self.args
  356. if argindex == 1:
  357. return (elliptic_e(m) + (m - n)*elliptic_k(m)/n +
  358. (n**2 - m)*elliptic_pi(n, m)/n)/(2*(m - n)*(n - 1))
  359. elif argindex == 2:
  360. return (elliptic_e(m)/(m - 1) + elliptic_pi(n, m))/(2*(n - m))
  361. raise ArgumentIndexError(self, argindex)
  362. def _eval_rewrite_as_Integral(self, *args, **kwargs):
  363. from sympy.integrals.integrals import Integral
  364. if len(self.args) == 2:
  365. n, m, z = self.args[0], self.args[1], pi/2
  366. else:
  367. n, z, m = self.args
  368. t = Dummy(uniquely_named_symbol('t', args).name)
  369. return Integral(1/((1 - n*sin(t)**2)*sqrt(1 - m*sin(t)**2)), (t, 0, z))