Coverage for /usr/lib/python3/dist-packages/sympy/functions/special/elliptic_integrals.py: 19%

1""" Elliptic Integrals. """ 

2 

3from sympy.core import S, pi, I, Rational 

4from sympy.core.function import Function, ArgumentIndexError 

5from sympy.core.symbol import Dummy 

6from sympy.functions.elementary.complexes import sign 

7from sympy.functions.elementary.hyperbolic import atanh 

8from sympy.functions.elementary.miscellaneous import sqrt 

9from sympy.functions.elementary.trigonometric import sin, tan 

10from sympy.functions.special.gamma_functions import gamma 

11from sympy.functions.special.hyper import hyper, meijerg 

12 

13class elliptic_k(Function): 

14 r""" 

15 The complete elliptic integral of the first kind, defined by 

16 

17 .. math:: K(m) = F\left(\tfrac{\pi}{2}\middle| m\right) 

18 

19 where $F\left(z\middle| m\right)$ is the Legendre incomplete 

20 elliptic integral of the first kind. 

21 

22 Explanation 

23 =========== 

24 

25 The function $K(m)$ is a single-valued function on the complex 

26 plane with branch cut along the interval $(1, \infty)$. 

27 

28 Note that our notation defines the incomplete elliptic integral 

29 in terms of the parameter $m$ instead of the elliptic modulus 

30 (eccentricity) $k$. 

31 In this case, the parameter $m$ is defined as $m=k^2$. 

32 

33 Examples 

34 ======== 

35 

36 >>> from sympy import elliptic_k, I 

37 >>> from sympy.abc import m 

38 >>> elliptic_k(0) 

39 pi/2 

40 >>> elliptic_k(1.0 + I) 

41 1.50923695405127 + 0.625146415202697*I 

42 >>> elliptic_k(m).series(n=3) 

43 pi/2 + pi*m/8 + 9*pi*m**2/128 + O(m**3) 

44 

45 See Also 

46 ======== 

47 

48 elliptic_f 

49 

50 References 

51 ========== 

52 

53 .. [1] https://en.wikipedia.org/wiki/Elliptic_integrals 

54 .. [2] https://functions.wolfram.com/EllipticIntegrals/EllipticK 

55 

56 """ 

57 

58 @classmethod 

59 def eval(cls, m): 

60 if m.is_zero: 

61 return pi*S.Half 

62 elif m is S.Half: 

63 return 8*pi**Rational(3, 2)/gamma(Rational(-1, 4))**2 

64 elif m is S.One: 

65 return S.ComplexInfinity 

66 elif m is S.NegativeOne: 

67 return gamma(Rational(1, 4))**2/(4*sqrt(2*pi)) 

68 elif m in (S.Infinity, S.NegativeInfinity, I*S.Infinity, 

69 I*S.NegativeInfinity, S.ComplexInfinity): 

70 return S.Zero 

71 

72 def fdiff(self, argindex=1): 

73 m = self.args[0] 

74 return (elliptic_e(m) - (1 - m)*elliptic_k(m))/(2*m*(1 - m)) 

75 

76 def _eval_conjugate(self): 

77 m = self.args[0] 

78 if (m.is_real and (m - 1).is_positive) is False: 

79 return self.func(m.conjugate()) 

80 

81 def _eval_nseries(self, x, n, logx, cdir=0): 

82 from sympy.simplify import hyperexpand 

83 return hyperexpand(self.rewrite(hyper)._eval_nseries(x, n=n, logx=logx)) 

84 

85 def _eval_rewrite_as_hyper(self, m, **kwargs): 

86 return pi*S.Half*hyper((S.Half, S.Half), (S.One,), m) 

87 

88 def _eval_rewrite_as_meijerg(self, m, **kwargs): 

89 return meijerg(((S.Half, S.Half), []), ((S.Zero,), (S.Zero,)), -m)/2 

90 

91 def _eval_is_zero(self): 

92 m = self.args[0] 

93 if m.is_infinite: 

94 return True 

95 

96 def _eval_rewrite_as_Integral(self, *args): 

97 from sympy.integrals.integrals import Integral 

98 t = Dummy('t') 

99 m = self.args[0] 

100 return Integral(1/sqrt(1 - m*sin(t)**2), (t, 0, pi/2)) 

101 

102 

103class elliptic_f(Function): 

104 r""" 

105 The Legendre incomplete elliptic integral of the first 

106 kind, defined by 

107 

108 .. math:: F\left(z\middle| m\right) = 

109 \int_0^z \frac{dt}{\sqrt{1 - m \sin^2 t}} 

110 

111 Explanation 

112 =========== 

113 

114 This function reduces to a complete elliptic integral of 

115 the first kind, $K(m)$, when $z = \pi/2$. 

116 

117 Note that our notation defines the incomplete elliptic integral 

118 in terms of the parameter $m$ instead of the elliptic modulus 

119 (eccentricity) $k$. 

120 In this case, the parameter $m$ is defined as $m=k^2$. 

121 

122 Examples 

123 ======== 

124 

125 >>> from sympy import elliptic_f, I 

126 >>> from sympy.abc import z, m 

127 >>> elliptic_f(z, m).series(z) 

128 z + z**5*(3*m**2/40 - m/30) + m*z**3/6 + O(z**6) 

129 >>> elliptic_f(3.0 + I/2, 1.0 + I) 

130 2.909449841483 + 1.74720545502474*I 

131 

132 See Also 

133 ======== 

134 

135 elliptic_k 

136 

137 References 

138 ========== 

139 

140 .. [1] https://en.wikipedia.org/wiki/Elliptic_integrals 

141 .. [2] https://functions.wolfram.com/EllipticIntegrals/EllipticF 

142 

143 """ 

144 

145 @classmethod 

146 def eval(cls, z, m): 

147 if z.is_zero: 

148 return S.Zero 

149 if m.is_zero: 

150 return z 

151 k = 2*z/pi 

152 if k.is_integer: 

153 return k*elliptic_k(m) 

154 elif m in (S.Infinity, S.NegativeInfinity): 

155 return S.Zero 

156 elif z.could_extract_minus_sign(): 

157 return -elliptic_f(-z, m) 

158 

159 def fdiff(self, argindex=1): 

160 z, m = self.args 

161 fm = sqrt(1 - m*sin(z)**2) 

162 if argindex == 1: 

163 return 1/fm 

164 elif argindex == 2: 

165 return (elliptic_e(z, m)/(2*m*(1 - m)) - elliptic_f(z, m)/(2*m) - 

166 sin(2*z)/(4*(1 - m)*fm)) 

167 raise ArgumentIndexError(self, argindex) 

168 

169 def _eval_conjugate(self): 

170 z, m = self.args 

171 if (m.is_real and (m - 1).is_positive) is False: 

172 return self.func(z.conjugate(), m.conjugate()) 

173 

174 def _eval_rewrite_as_Integral(self, *args): 

175 from sympy.integrals.integrals import Integral 

176 t = Dummy('t') 

177 z, m = self.args[0], self.args[1] 

178 return Integral(1/(sqrt(1 - m*sin(t)**2)), (t, 0, z)) 

179 

180 def _eval_is_zero(self): 

181 z, m = self.args 

182 if z.is_zero: 

183 return True 

184 if m.is_extended_real and m.is_infinite: 

185 return True 

186 

187 

188class elliptic_e(Function): 

189 r""" 

190 Called with two arguments $z$ and $m$, evaluates the 

191 incomplete elliptic integral of the second kind, defined by 

192 

193 .. math:: E\left(z\middle| m\right) = \int_0^z \sqrt{1 - m \sin^2 t} dt 

194 

195 Called with a single argument $m$, evaluates the Legendre complete 

196 elliptic integral of the second kind 

197 

198 .. math:: E(m) = E\left(\tfrac{\pi}{2}\middle| m\right) 

199 

200 Explanation 

201 =========== 

202 

203 The function $E(m)$ is a single-valued function on the complex 

204 plane with branch cut along the interval $(1, \infty)$. 

205 

206 Note that our notation defines the incomplete elliptic integral 

207 in terms of the parameter $m$ instead of the elliptic modulus 

208 (eccentricity) $k$. 

209 In this case, the parameter $m$ is defined as $m=k^2$. 

210 

211 Examples 

212 ======== 

213 

214 >>> from sympy import elliptic_e, I 

215 >>> from sympy.abc import z, m 

216 >>> elliptic_e(z, m).series(z) 

217 z + z**5*(-m**2/40 + m/30) - m*z**3/6 + O(z**6) 

218 >>> elliptic_e(m).series(n=4) 

219 pi/2 - pi*m/8 - 3*pi*m**2/128 - 5*pi*m**3/512 + O(m**4) 

220 >>> elliptic_e(1 + I, 2 - I/2).n() 

221 1.55203744279187 + 0.290764986058437*I 

222 >>> elliptic_e(0) 

223 pi/2 

224 >>> elliptic_e(2.0 - I) 

225 0.991052601328069 + 0.81879421395609*I 

226 

227 References 

228 ========== 

229 

230 .. [1] https://en.wikipedia.org/wiki/Elliptic_integrals 

231 .. [2] https://functions.wolfram.com/EllipticIntegrals/EllipticE2 

232 .. [3] https://functions.wolfram.com/EllipticIntegrals/EllipticE 

233 

234 """ 

235 

236 @classmethod 

237 def eval(cls, m, z=None): 

238 if z is not None: 

239 z, m = m, z 

240 k = 2*z/pi 

241 if m.is_zero: 

242 return z 

243 if z.is_zero: 

244 return S.Zero 

245 elif k.is_integer: 

246 return k*elliptic_e(m) 

247 elif m in (S.Infinity, S.NegativeInfinity): 

248 return S.ComplexInfinity 

249 elif z.could_extract_minus_sign(): 

250 return -elliptic_e(-z, m) 

251 else: 

252 if m.is_zero: 

253 return pi/2 

254 elif m is S.One: 

255 return S.One 

256 elif m is S.Infinity: 

257 return I*S.Infinity 

258 elif m is S.NegativeInfinity: 

259 return S.Infinity 

260 elif m is S.ComplexInfinity: 

261 return S.ComplexInfinity 

262 

263 def fdiff(self, argindex=1): 

264 if len(self.args) == 2: 

265 z, m = self.args 

266 if argindex == 1: 

267 return sqrt(1 - m*sin(z)**2) 

268 elif argindex == 2: 

269 return (elliptic_e(z, m) - elliptic_f(z, m))/(2*m) 

270 else: 

271 m = self.args[0] 

272 if argindex == 1: 

273 return (elliptic_e(m) - elliptic_k(m))/(2*m) 

274 raise ArgumentIndexError(self, argindex) 

275 

276 def _eval_conjugate(self): 

277 if len(self.args) == 2: 

278 z, m = self.args 

279 if (m.is_real and (m - 1).is_positive) is False: 

280 return self.func(z.conjugate(), m.conjugate()) 

281 else: 

282 m = self.args[0] 

283 if (m.is_real and (m - 1).is_positive) is False: 

284 return self.func(m.conjugate()) 

285 

286 def _eval_nseries(self, x, n, logx, cdir=0): 

287 from sympy.simplify import hyperexpand 

288 if len(self.args) == 1: 

289 return hyperexpand(self.rewrite(hyper)._eval_nseries(x, n=n, logx=logx)) 

290 return super()._eval_nseries(x, n=n, logx=logx) 

291 

292 def _eval_rewrite_as_hyper(self, *args, **kwargs): 

293 if len(args) == 1: 

294 m = args[0] 

295 return (pi/2)*hyper((Rational(-1, 2), S.Half), (S.One,), m) 

296 

297 def _eval_rewrite_as_meijerg(self, *args, **kwargs): 

298 if len(args) == 1: 

299 m = args[0] 

300 return -meijerg(((S.Half, Rational(3, 2)), []), \ 

301 ((S.Zero,), (S.Zero,)), -m)/4 

302 

303 def _eval_rewrite_as_Integral(self, *args): 

304 from sympy.integrals.integrals import Integral 

305 z, m = (pi/2, self.args[0]) if len(self.args) == 1 else self.args 

306 t = Dummy('t') 

307 return Integral(sqrt(1 - m*sin(t)**2), (t, 0, z)) 

308 

309 

310class elliptic_pi(Function): 

311 r""" 

312 Called with three arguments $n$, $z$ and $m$, evaluates the 

313 Legendre incomplete elliptic integral of the third kind, defined by 

314 

315 .. math:: \Pi\left(n; z\middle| m\right) = \int_0^z \frac{dt} 

316 {\left(1 - n \sin^2 t\right) \sqrt{1 - m \sin^2 t}} 

317 

318 Called with two arguments $n$ and $m$, evaluates the complete 

319 elliptic integral of the third kind: 

320 

321 .. math:: \Pi\left(n\middle| m\right) = 

322 \Pi\left(n; \tfrac{\pi}{2}\middle| m\right) 

323 

324 Explanation 

325 =========== 

326 

327 Note that our notation defines the incomplete elliptic integral 

328 in terms of the parameter $m$ instead of the elliptic modulus 

329 (eccentricity) $k$. 

330 In this case, the parameter $m$ is defined as $m=k^2$. 

331 

332 Examples 

333 ======== 

334 

335 >>> from sympy import elliptic_pi, I 

336 >>> from sympy.abc import z, n, m 

337 >>> elliptic_pi(n, z, m).series(z, n=4) 

338 z + z**3*(m/6 + n/3) + O(z**4) 

339 >>> elliptic_pi(0.5 + I, 1.0 - I, 1.2) 

340 2.50232379629182 - 0.760939574180767*I 

341 >>> elliptic_pi(0, 0) 

342 pi/2 

343 >>> elliptic_pi(1.0 - I/3, 2.0 + I) 

344 3.29136443417283 + 0.32555634906645*I 

345 

346 References 

347 ========== 

348 

349 .. [1] https://en.wikipedia.org/wiki/Elliptic_integrals 

350 .. [2] https://functions.wolfram.com/EllipticIntegrals/EllipticPi3 

351 .. [3] https://functions.wolfram.com/EllipticIntegrals/EllipticPi 

352 

353 """ 

354 

355 @classmethod 

356 def eval(cls, n, m, z=None): 

357 if z is not None: 

358 n, z, m = n, m, z 

359 if n.is_zero: 

360 return elliptic_f(z, m) 

361 elif n is S.One: 

362 return (elliptic_f(z, m) + 

363 (sqrt(1 - m*sin(z)**2)*tan(z) - 

364 elliptic_e(z, m))/(1 - m)) 

365 k = 2*z/pi 

366 if k.is_integer: 

367 return k*elliptic_pi(n, m) 

368 elif m.is_zero: 

369 return atanh(sqrt(n - 1)*tan(z))/sqrt(n - 1) 

370 elif n == m: 

371 return (elliptic_f(z, n) - elliptic_pi(1, z, n) + 

372 tan(z)/sqrt(1 - n*sin(z)**2)) 

373 elif n in (S.Infinity, S.NegativeInfinity): 

374 return S.Zero 

375 elif m in (S.Infinity, S.NegativeInfinity): 

376 return S.Zero 

377 elif z.could_extract_minus_sign(): 

378 return -elliptic_pi(n, -z, m) 

379 if n.is_zero: 

380 return elliptic_f(z, m) 

381 if m.is_extended_real and m.is_infinite or \ 

382 n.is_extended_real and n.is_infinite: 

383 return S.Zero 

384 else: 

385 if n.is_zero: 

386 return elliptic_k(m) 

387 elif n is S.One: 

388 return S.ComplexInfinity 

389 elif m.is_zero: 

390 return pi/(2*sqrt(1 - n)) 

391 elif m == S.One: 

392 return S.NegativeInfinity/sign(n - 1) 

393 elif n == m: 

394 return elliptic_e(n)/(1 - n) 

395 elif n in (S.Infinity, S.NegativeInfinity): 

396 return S.Zero 

397 elif m in (S.Infinity, S.NegativeInfinity): 

398 return S.Zero 

399 if n.is_zero: 

400 return elliptic_k(m) 

401 if m.is_extended_real and m.is_infinite or \ 

402 n.is_extended_real and n.is_infinite: 

403 return S.Zero 

404 

405 def _eval_conjugate(self): 

406 if len(self.args) == 3: 

407 n, z, m = self.args 

408 if (n.is_real and (n - 1).is_positive) is False and \ 

409 (m.is_real and (m - 1).is_positive) is False: 

410 return self.func(n.conjugate(), z.conjugate(), m.conjugate()) 

411 else: 

412 n, m = self.args 

413 return self.func(n.conjugate(), m.conjugate()) 

414 

415 def fdiff(self, argindex=1): 

416 if len(self.args) == 3: 

417 n, z, m = self.args 

418 fm, fn = sqrt(1 - m*sin(z)**2), 1 - n*sin(z)**2 

419 if argindex == 1: 

420 return (elliptic_e(z, m) + (m - n)*elliptic_f(z, m)/n + 

421 (n**2 - m)*elliptic_pi(n, z, m)/n - 

422 n*fm*sin(2*z)/(2*fn))/(2*(m - n)*(n - 1)) 

423 elif argindex == 2: 

424 return 1/(fm*fn) 

425 elif argindex == 3: 

426 return (elliptic_e(z, m)/(m - 1) + 

427 elliptic_pi(n, z, m) - 

428 m*sin(2*z)/(2*(m - 1)*fm))/(2*(n - m)) 

429 else: 

430 n, m = self.args 

431 if argindex == 1: 

432 return (elliptic_e(m) + (m - n)*elliptic_k(m)/n + 

433 (n**2 - m)*elliptic_pi(n, m)/n)/(2*(m - n)*(n - 1)) 

434 elif argindex == 2: 

435 return (elliptic_e(m)/(m - 1) + elliptic_pi(n, m))/(2*(n - m)) 

436 raise ArgumentIndexError(self, argindex) 

437 

438 def _eval_rewrite_as_Integral(self, *args): 

439 from sympy.integrals.integrals import Integral 

440 if len(self.args) == 2: 

441 n, m, z = self.args[0], self.args[1], pi/2 

442 else: 

443 n, z, m = self.args 

444 t = Dummy('t') 

445 return Integral(1/((1 - n*sin(t)**2)*sqrt(1 - m*sin(t)**2)), (t, 0, z))