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

1from functools import wraps

3from sympy.core import S

4from sympy.core.add import Add

5from sympy.core.cache import cacheit

6from sympy.core.expr import Expr

7from sympy.core.function import Function, ArgumentIndexError, _mexpand

8from sympy.core.logic import fuzzy_or, fuzzy_not

9from sympy.core.numbers import Rational, pi, I

10from sympy.core.power import Pow

11from sympy.core.symbol import Dummy, Wild

12from sympy.core.sympify import sympify

13from sympy.functions.combinatorial.factorials import factorial

14from sympy.functions.elementary.trigonometric import sin, cos, csc, cot

15from sympy.functions.elementary.integers import ceiling

16from sympy.functions.elementary.exponential import exp, log

17from sympy.functions.elementary.miscellaneous import cbrt, sqrt, root

18from sympy.functions.elementary.complexes import (Abs, re, im, polar_lift, unpolarify)

19from sympy.functions.special.gamma_functions import gamma, digamma, uppergamma

20from sympy.functions.special.hyper import hyper

21from sympy.polys.orthopolys import spherical_bessel_fn

23from mpmath import mp, workprec

25# TODO

26# o Scorer functions G1 and G2

27# o Asymptotic expansions

28# These are possible, e.g. for fixed order, but since the bessel type

29# functions are oscillatory they are not actually tractable at

30# infinity, so this is not particularly useful right now.

31# o Nicer series expansions.

32# o More rewriting.

33# o Add solvers to ode.py (or rather add solvers for the hypergeometric equation).

36class BesselBase(Function):

37 """

38 Abstract base class for Bessel-type functions.

40 This class is meant to reduce code duplication.

41 All Bessel-type functions can 1) be differentiated, with the derivatives

42 expressed in terms of similar functions, and 2) be rewritten in terms

43 of other Bessel-type functions.

45 Here, Bessel-type functions are assumed to have one complex parameter.

47 To use this base class, define class attributes ``_a`` and ``_b`` such that

48 ``2*F_n' = -_a*F_{n+1} + b*F_{n-1}``.

50 """

52 @property

53 def order(self):

54 """ The order of the Bessel-type function. """

55 return self.args[0]

57 @property

58 def argument(self):

59 """ The argument of the Bessel-type function. """

60 return self.args[1]

62 @classmethod

63 def eval(cls, nu, z):

64 return

66 def fdiff(self, argindex=2):

67 if argindex != 2:

68 raise ArgumentIndexError(self, argindex)

69 return (self._b/2 * self.__class__(self.order - 1, self.argument) -

70 self._a/2 * self.__class__(self.order + 1, self.argument))

72 def _eval_conjugate(self):

73 z = self.argument

74 if z.is_extended_negative is False:

75 return self.__class__(self.order.conjugate(), z.conjugate())

77 def _eval_is_meromorphic(self, x, a):

78 nu, z = self.order, self.argument

80 if nu.has(x):

81 return False

82 if not z._eval_is_meromorphic(x, a):

83 return None

84 z0 = z.subs(x, a)

85 if nu.is_integer:

86 if isinstance(self, (besselj, besseli, hn1, hn2, jn, yn)) or not nu.is_zero:

87 return fuzzy_not(z0.is_infinite)

88 return fuzzy_not(fuzzy_or([z0.is_zero, z0.is_infinite]))

90 def _eval_expand_func(self, **hints):

91 nu, z, f = self.order, self.argument, self.__class__

92 if nu.is_real:

93 if (nu - 1).is_positive:

94 return (-self._a*self._b*f(nu - 2, z)._eval_expand_func() +

95 2*self._a*(nu - 1)*f(nu - 1, z)._eval_expand_func()/z)

96 elif (nu + 1).is_negative:

97 return (2*self._b*(nu + 1)*f(nu + 1, z)._eval_expand_func()/z -

98 self._a*self._b*f(nu + 2, z)._eval_expand_func())

99 return self

100

101 def _eval_simplify(self, **kwargs):

102 from sympy.simplify.simplify import besselsimp

103 return besselsimp(self)

104

105

106class besselj(BesselBase):

107 r"""

108 Bessel function of the first kind.

109

110 Explanation

111 ===========

112

113 The Bessel $J$ function of order $\nu$ is defined to be the function

114 satisfying Bessel's differential equation

115

116 .. math ::

117 z^2 \frac{\mathrm{d}^2 w}{\mathrm{d}z^2}

118 + z \frac{\mathrm{d}w}{\mathrm{d}z} + (z^2 - \nu^2) w = 0,

119

120 with Laurent expansion

121

122 .. math ::

123 J_\nu(z) = z^\nu \left(\frac{1}{\Gamma(\nu + 1) 2^\nu} + O(z^2) \right),

124

125 if $\nu$ is not a negative integer. If $\nu=-n \in \mathbb{Z}_{<0}$

126 *is* a negative integer, then the definition is

127

128 .. math ::

129 J_{-n}(z) = (-1)^n J_n(z).

130

131 Examples

132 ========

133

134 Create a Bessel function object:

135

136 >>> from sympy import besselj, jn

137 >>> from sympy.abc import z, n

138 >>> b = besselj(n, z)

139

140 Differentiate it:

141

142 >>> b.diff(z)

143 besselj(n - 1, z)/2 - besselj(n + 1, z)/2

144

145 Rewrite in terms of spherical Bessel functions:

146

147 >>> b.rewrite(jn)

148 sqrt(2)*sqrt(z)*jn(n - 1/2, z)/sqrt(pi)

149

150 Access the parameter and argument:

151

152 >>> b.order

153 n

154 >>> b.argument

155 z

156

157 See Also

158 ========

159

160 bessely, besseli, besselk

161

162 References

163 ==========

164

165 .. [1] Abramowitz, Milton; Stegun, Irene A., eds. (1965), "Chapter 9",

166 Handbook of Mathematical Functions with Formulas, Graphs, and

167 Mathematical Tables

168 .. [2] Luke, Y. L. (1969), The Special Functions and Their

169 Approximations, Volume 1

170 .. [3] https://en.wikipedia.org/wiki/Bessel_function

171 .. [4] https://functions.wolfram.com/Bessel-TypeFunctions/BesselJ/

172

173 """

174

175 _a = S.One

176 _b = S.One

177

178 @classmethod

179 def eval(cls, nu, z):

180 if z.is_zero:

181 if nu.is_zero:

182 return S.One

183 elif (nu.is_integer and nu.is_zero is False) or re(nu).is_positive:

184 return S.Zero

185 elif re(nu).is_negative and not (nu.is_integer is True):

186 return S.ComplexInfinity

187 elif nu.is_imaginary:

188 return S.NaN

189 if z in (S.Infinity, S.NegativeInfinity):

190 return S.Zero

191

192 if z.could_extract_minus_sign():

193 return (z)**nu*(-z)**(-nu)*besselj(nu, -z)

194 if nu.is_integer:

195 if nu.could_extract_minus_sign():

196 return S.NegativeOne**(-nu)*besselj(-nu, z)

197 newz = z.extract_multiplicatively(I)

198 if newz: # NOTE we don't want to change the function if z==0

199 return I**(nu)*besseli(nu, newz)

200

201 # branch handling:

202 if nu.is_integer:

203 newz = unpolarify(z)

204 if newz != z:

205 return besselj(nu, newz)

206 else:

207 newz, n = z.extract_branch_factor()

208 if n != 0:

209 return exp(2*n*pi*nu*I)*besselj(nu, newz)

210 nnu = unpolarify(nu)

211 if nu != nnu:

212 return besselj(nnu, z)

213

214 def _eval_rewrite_as_besseli(self, nu, z, **kwargs):

215 return exp(I*pi*nu/2)*besseli(nu, polar_lift(-I)*z)

216

217 def _eval_rewrite_as_bessely(self, nu, z, **kwargs):

218 if nu.is_integer is False:

219 return csc(pi*nu)*bessely(-nu, z) - cot(pi*nu)*bessely(nu, z)

220

221 def _eval_rewrite_as_jn(self, nu, z, **kwargs):

222 return sqrt(2*z/pi)*jn(nu - S.Half, self.argument)

223

224 def _eval_as_leading_term(self, x, logx=None, cdir=0):

225 nu, z = self.args

226 try:

227 arg = z.as_leading_term(x)

228 except NotImplementedError:

229 return self

230 c, e = arg.as_coeff_exponent(x)

231

232 if e.is_positive:

233 return arg**nu/(2**nu*gamma(nu + 1))

234 elif e.is_negative:

235 cdir = 1 if cdir == 0 else cdir

236 sign = c*cdir**e

237 if not sign.is_negative:

238 # Refer Abramowitz and Stegun 1965, p. 364 for more information on

239 # asymptotic approximation of besselj function.

240 return sqrt(2)*cos(z - pi*(2*nu + 1)/4)/sqrt(pi*z)

241 return self

242

243 return super(besselj, self)._eval_as_leading_term(x, logx, cdir)

244

245 def _eval_is_extended_real(self):

246 nu, z = self.args

247 if nu.is_integer and z.is_extended_real:

248 return True

249

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

251 # Refer https://functions.wolfram.com/Bessel-TypeFunctions/BesselJ/06/01/04/01/01/0003/

252 # for more information on nseries expansion of besselj function.

253 from sympy.series.order import Order

254 nu, z = self.args

255

256 # In case of powers less than 1, number of terms need to be computed

257 # separately to avoid repeated callings of _eval_nseries with wrong n

258 try:

259 _, exp = z.leadterm(x)

260 except (ValueError, NotImplementedError):

261 return self

262

263 if exp.is_positive:

264 newn = ceiling(n/exp)

265 o = Order(x**n, x)

266 r = (z/2)._eval_nseries(x, n, logx, cdir).removeO()

267 if r is S.Zero:

268 return o

269 t = (_mexpand(r**2) + o).removeO()

270

271 term = r**nu/gamma(nu + 1)

272 s = [term]

273 for k in range(1, (newn + 1)//2):

274 term *= -t/(k*(nu + k))

275 term = (_mexpand(term) + o).removeO()

276 s.append(term)

277 return Add(*s) + o

278

279 return super(besselj, self)._eval_nseries(x, n, logx, cdir)

280

281

282class bessely(BesselBase):

283 r"""

284 Bessel function of the second kind.

285

286 Explanation

287 ===========

288

289 The Bessel $Y$ function of order $\nu$ is defined as

290

291 .. math ::

292 Y_\nu(z) = \lim_{\mu \to \nu} \frac{J_\mu(z) \cos(\pi \mu)

293 - J_{-\mu}(z)}{\sin(\pi \mu)},

294

295 where $J_\mu(z)$ is the Bessel function of the first kind.

296

297 It is a solution to Bessel's equation, and linearly independent from

298 $J_\nu$.

299

300 Examples

301 ========

302

303 >>> from sympy import bessely, yn

304 >>> from sympy.abc import z, n

305 >>> b = bessely(n, z)

306 >>> b.diff(z)

307 bessely(n - 1, z)/2 - bessely(n + 1, z)/2

308 >>> b.rewrite(yn)

309 sqrt(2)*sqrt(z)*yn(n - 1/2, z)/sqrt(pi)

310

311 See Also

312 ========

313

314 besselj, besseli, besselk

315

316 References

317 ==========

318

319 .. [1] https://functions.wolfram.com/Bessel-TypeFunctions/BesselY/

320

321 """

322

323 _a = S.One

324 _b = S.One

325

326 @classmethod

327 def eval(cls, nu, z):

328 if z.is_zero:

329 if nu.is_zero:

330 return S.NegativeInfinity

331 elif re(nu).is_zero is False:

332 return S.ComplexInfinity

333 elif re(nu).is_zero:

334 return S.NaN

335 if z in (S.Infinity, S.NegativeInfinity):

336 return S.Zero

337 if z == I*S.Infinity:

338 return exp(I*pi*(nu + 1)/2) * S.Infinity

339 if z == I*S.NegativeInfinity:

340 return exp(-I*pi*(nu + 1)/2) * S.Infinity

341

342 if nu.is_integer:

343 if nu.could_extract_minus_sign():

344 return S.NegativeOne**(-nu)*bessely(-nu, z)

345

346 def _eval_rewrite_as_besselj(self, nu, z, **kwargs):

347 if nu.is_integer is False:

348 return csc(pi*nu)*(cos(pi*nu)*besselj(nu, z) - besselj(-nu, z))

349

350 def _eval_rewrite_as_besseli(self, nu, z, **kwargs):

351 aj = self._eval_rewrite_as_besselj(*self.args)

352 if aj:

353 return aj.rewrite(besseli)

354

355 def _eval_rewrite_as_yn(self, nu, z, **kwargs):

356 return sqrt(2*z/pi) * yn(nu - S.Half, self.argument)

357

358 def _eval_as_leading_term(self, x, logx=None, cdir=0):

359 nu, z = self.args

360 try:

361 arg = z.as_leading_term(x)

362 except NotImplementedError:

363 return self

364 c, e = arg.as_coeff_exponent(x)

365

366 if e.is_positive:

367 term_one = ((2/pi)*log(z/2)*besselj(nu, z))

368 term_two = -(z/2)**(-nu)*factorial(nu - 1)/pi if (nu).is_positive else S.Zero

369 term_three = -(z/2)**nu/(pi*factorial(nu))*(digamma(nu + 1) - S.EulerGamma)

370 arg = Add(*[term_one, term_two, term_three]).as_leading_term(x, logx=logx)

371 return arg

372 elif e.is_negative:

373 cdir = 1 if cdir == 0 else cdir

374 sign = c*cdir**e

375 if not sign.is_negative:

376 # Refer Abramowitz and Stegun 1965, p. 364 for more information on

377 # asymptotic approximation of bessely function.

378 return sqrt(2)*(-sin(pi*nu/2 - z + pi/4) + 3*cos(pi*nu/2 - z + pi/4)/(8*z))*sqrt(1/z)/sqrt(pi)

379 return self

380

381 return super(bessely, self)._eval_as_leading_term(x, logx, cdir)

382

383 def _eval_is_extended_real(self):

384 nu, z = self.args

385 if nu.is_integer and z.is_positive:

386 return True

387

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

389 # Refer https://functions.wolfram.com/Bessel-TypeFunctions/BesselY/06/01/04/01/02/0008/

390 # for more information on nseries expansion of bessely function.

391 from sympy.series.order import Order

392 nu, z = self.args

393

394 # In case of powers less than 1, number of terms need to be computed

395 # separately to avoid repeated callings of _eval_nseries with wrong n

396 try:

397 _, exp = z.leadterm(x)

398 except (ValueError, NotImplementedError):

399 return self

400

401 if exp.is_positive and nu.is_integer:

402 newn = ceiling(n/exp)

403 bn = besselj(nu, z)

404 a = ((2/pi)*log(z/2)*bn)._eval_nseries(x, n, logx, cdir)

405

406 b, c = [], []

407 o = Order(x**n, x)

408 r = (z/2)._eval_nseries(x, n, logx, cdir).removeO()

409 if r is S.Zero:

410 return o

411 t = (_mexpand(r**2) + o).removeO()

412

413 if nu > S.Zero:

414 term = r**(-nu)*factorial(nu - 1)/pi

415 b.append(term)

416 for k in range(1, nu):

417 denom = (nu - k)*k

418 if denom == S.Zero:

419 term *= t/k

420 else:

421 term *= t/denom

422 term = (_mexpand(term) + o).removeO()

423 b.append(term)

424

425 p = r**nu/(pi*factorial(nu))

426 term = p*(digamma(nu + 1) - S.EulerGamma)

427 c.append(term)

428 for k in range(1, (newn + 1)//2):

429 p *= -t/(k*(k + nu))

430 p = (_mexpand(p) + o).removeO()

431 term = p*(digamma(k + nu + 1) + digamma(k + 1))

432 c.append(term)

433 return a - Add(*b) - Add(*c) # Order term comes from a

434

435 return super(bessely, self)._eval_nseries(x, n, logx, cdir)

436

437

438class besseli(BesselBase):

439 r"""

440 Modified Bessel function of the first kind.

441

442 Explanation

443 ===========

444

445 The Bessel $I$ function is a solution to the modified Bessel equation

446

447 .. math ::

448 z^2 \frac{\mathrm{d}^2 w}{\mathrm{d}z^2}

449 + z \frac{\mathrm{d}w}{\mathrm{d}z} + (z^2 + \nu^2)^2 w = 0.

450

451 It can be defined as

452

453 .. math ::

454 I_\nu(z) = i^{-\nu} J_\nu(iz),

455

456 where $J_\nu(z)$ is the Bessel function of the first kind.

457

458 Examples

459 ========

460

461 >>> from sympy import besseli

462 >>> from sympy.abc import z, n

463 >>> besseli(n, z).diff(z)

464 besseli(n - 1, z)/2 + besseli(n + 1, z)/2

465

466 See Also

467 ========

468

469 besselj, bessely, besselk

470

471 References

472 ==========

473

474 .. [1] https://functions.wolfram.com/Bessel-TypeFunctions/BesselI/

475

476 """

477

478 _a = -S.One

479 _b = S.One

480

481 @classmethod

482 def eval(cls, nu, z):

483 if z.is_zero:

484 if nu.is_zero:

485 return S.One

486 elif (nu.is_integer and nu.is_zero is False) or re(nu).is_positive:

487 return S.Zero

488 elif re(nu).is_negative and not (nu.is_integer is True):

489 return S.ComplexInfinity

490 elif nu.is_imaginary:

491 return S.NaN

492 if im(z) in (S.Infinity, S.NegativeInfinity):

493 return S.Zero

494 if z is S.Infinity:

495 return S.Infinity

496 if z is S.NegativeInfinity:

497 return (-1)**nu*S.Infinity

498

499 if z.could_extract_minus_sign():

500 return (z)**nu*(-z)**(-nu)*besseli(nu, -z)

501 if nu.is_integer:

502 if nu.could_extract_minus_sign():

503 return besseli(-nu, z)

504 newz = z.extract_multiplicatively(I)

505 if newz: # NOTE we don't want to change the function if z==0

506 return I**(-nu)*besselj(nu, -newz)

507

508 # branch handling:

509 if nu.is_integer:

510 newz = unpolarify(z)

511 if newz != z:

512 return besseli(nu, newz)

513 else:

514 newz, n = z.extract_branch_factor()

515 if n != 0:

516 return exp(2*n*pi*nu*I)*besseli(nu, newz)

517 nnu = unpolarify(nu)

518 if nu != nnu:

519 return besseli(nnu, z)

520

521 def _eval_rewrite_as_besselj(self, nu, z, **kwargs):

522 return exp(-I*pi*nu/2)*besselj(nu, polar_lift(I)*z)

523

524 def _eval_rewrite_as_bessely(self, nu, z, **kwargs):

525 aj = self._eval_rewrite_as_besselj(*self.args)

526 if aj:

527 return aj.rewrite(bessely)

528

529 def _eval_rewrite_as_jn(self, nu, z, **kwargs):

530 return self._eval_rewrite_as_besselj(*self.args).rewrite(jn)

531

532 def _eval_is_extended_real(self):

533 nu, z = self.args

534 if nu.is_integer and z.is_extended_real:

535 return True

536

537 def _eval_as_leading_term(self, x, logx=None, cdir=0):

538 nu, z = self.args

539 try:

540 arg = z.as_leading_term(x)

541 except NotImplementedError:

542 return self

543 c, e = arg.as_coeff_exponent(x)

544

545 if e.is_positive:

546 return arg**nu/(2**nu*gamma(nu + 1))

547 elif e.is_negative:

548 cdir = 1 if cdir == 0 else cdir

549 sign = c*cdir**e

550 if not sign.is_negative:

551 # Refer Abramowitz and Stegun 1965, p. 377 for more information on

552 # asymptotic approximation of besseli function.

553 return exp(z)/sqrt(2*pi*z)

554 return self

555

556 return super(besseli, self)._eval_as_leading_term(x, logx, cdir)

557

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

559 # Refer https://functions.wolfram.com/Bessel-TypeFunctions/BesselI/06/01/04/01/01/0003/

560 # for more information on nseries expansion of besseli function.

561 from sympy.series.order import Order

562 nu, z = self.args

563

564 # In case of powers less than 1, number of terms need to be computed

565 # separately to avoid repeated callings of _eval_nseries with wrong n

566 try:

567 _, exp = z.leadterm(x)

568 except (ValueError, NotImplementedError):

569 return self

570

571 if exp.is_positive:

572 newn = ceiling(n/exp)

573 o = Order(x**n, x)

574 r = (z/2)._eval_nseries(x, n, logx, cdir).removeO()

575 if r is S.Zero:

576 return o

577 t = (_mexpand(r**2) + o).removeO()

578

579 term = r**nu/gamma(nu + 1)

580 s = [term]

581 for k in range(1, (newn + 1)//2):

582 term *= t/(k*(nu + k))

583 term = (_mexpand(term) + o).removeO()

584 s.append(term)

585 return Add(*s) + o

586

587 return super(besseli, self)._eval_nseries(x, n, logx, cdir)

588

589

590class besselk(BesselBase):

591 r"""

592 Modified Bessel function of the second kind.

593

594 Explanation

595 ===========

596

597 The Bessel $K$ function of order $\nu$ is defined as

598

599 .. math ::

600 K_\nu(z) = \lim_{\mu \to \nu} \frac{\pi}{2}

601 \frac{I_{-\mu}(z) -I_\mu(z)}{\sin(\pi \mu)},

602

603 where $I_\mu(z)$ is the modified Bessel function of the first kind.

604

605 It is a solution of the modified Bessel equation, and linearly independent

606 from $Y_\nu$.

607

608 Examples

609 ========

610

611 >>> from sympy import besselk

612 >>> from sympy.abc import z, n

613 >>> besselk(n, z).diff(z)

614 -besselk(n - 1, z)/2 - besselk(n + 1, z)/2

615

616 See Also

617 ========

618

619 besselj, besseli, bessely

620

621 References

622 ==========

623

624 .. [1] https://functions.wolfram.com/Bessel-TypeFunctions/BesselK/

625

626 """

627

628 _a = S.One

629 _b = -S.One

630

631 @classmethod

632 def eval(cls, nu, z):

633 if z.is_zero:

634 if nu.is_zero:

635 return S.Infinity

636 elif re(nu).is_zero is False:

637 return S.ComplexInfinity

638 elif re(nu).is_zero:

639 return S.NaN

640 if z in (S.Infinity, I*S.Infinity, I*S.NegativeInfinity):

641 return S.Zero

642

643 if nu.is_integer:

644 if nu.could_extract_minus_sign():

645 return besselk(-nu, z)

646

647 def _eval_rewrite_as_besseli(self, nu, z, **kwargs):

648 if nu.is_integer is False:

649 return pi*csc(pi*nu)*(besseli(-nu, z) - besseli(nu, z))/2

650

651 def _eval_rewrite_as_besselj(self, nu, z, **kwargs):

652 ai = self._eval_rewrite_as_besseli(*self.args)

653 if ai:

654 return ai.rewrite(besselj)

655

656 def _eval_rewrite_as_bessely(self, nu, z, **kwargs):

657 aj = self._eval_rewrite_as_besselj(*self.args)

658 if aj:

659 return aj.rewrite(bessely)

660

661 def _eval_rewrite_as_yn(self, nu, z, **kwargs):

662 ay = self._eval_rewrite_as_bessely(*self.args)

663 if ay:

664 return ay.rewrite(yn)

665

666 def _eval_is_extended_real(self):

667 nu, z = self.args

668 if nu.is_integer and z.is_positive:

669 return True

670

671 def _eval_as_leading_term(self, x, logx=None, cdir=0):

672 nu, z = self.args

673 try:

674 arg = z.as_leading_term(x)

675 except NotImplementedError:

676 return self

677 _, e = arg.as_coeff_exponent(x)

678

679 if e.is_positive:

680 term_one = ((-1)**(nu -1)*log(z/2)*besseli(nu, z))

681 term_two = (z/2)**(-nu)*factorial(nu - 1)/2 if (nu).is_positive else S.Zero

682 term_three = (-1)**nu*(z/2)**nu/(2*factorial(nu))*(digamma(nu + 1) - S.EulerGamma)

683 arg = Add(*[term_one, term_two, term_three]).as_leading_term(x, logx=logx)

684 return arg

685 elif e.is_negative:

686 # Refer Abramowitz and Stegun 1965, p. 378 for more information on

687 # asymptotic approximation of besselk function.

688 return sqrt(pi)*exp(-z)/sqrt(2*z)

689

690 return super(besselk, self)._eval_as_leading_term(x, logx, cdir)

691

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

693 # Refer https://functions.wolfram.com/Bessel-TypeFunctions/BesselK/06/01/04/01/02/0008/

694 # for more information on nseries expansion of besselk function.

695 from sympy.series.order import Order

696 nu, z = self.args

697

698 # In case of powers less than 1, number of terms need to be computed

699 # separately to avoid repeated callings of _eval_nseries with wrong n

700 try:

701 _, exp = z.leadterm(x)

702 except (ValueError, NotImplementedError):

703 return self

704

705 if exp.is_positive and nu.is_integer:

706 newn = ceiling(n/exp)

707 bn = besseli(nu, z)

708 a = ((-1)**(nu - 1)*log(z/2)*bn)._eval_nseries(x, n, logx, cdir)

709

710 b, c = [], []

711 o = Order(x**n, x)

712 r = (z/2)._eval_nseries(x, n, logx, cdir).removeO()

713 if r is S.Zero:

714 return o

715 t = (_mexpand(r**2) + o).removeO()

716

717 if nu > S.Zero:

718 term = r**(-nu)*factorial(nu - 1)/2

719 b.append(term)

720 for k in range(1, nu):

721 denom = (k - nu)*k

722 if denom == S.Zero:

723 term *= t/k

724 else:

725 term *= t/denom

726 term = (_mexpand(term) + o).removeO()

727 b.append(term)

728

729 p = r**nu*(-1)**nu/(2*factorial(nu))

730 term = p*(digamma(nu + 1) - S.EulerGamma)

731 c.append(term)

732 for k in range(1, (newn + 1)//2):

733 p *= t/(k*(k + nu))

734 p = (_mexpand(p) + o).removeO()

735 term = p*(digamma(k + nu + 1) + digamma(k + 1))

736 c.append(term)

737 return a + Add(*b) + Add(*c) # Order term comes from a

738

739 return super(besselk, self)._eval_nseries(x, n, logx, cdir)

740

741

742class hankel1(BesselBase):

743 r"""

744 Hankel function of the first kind.

745

746 Explanation

747 ===========

748

749 This function is defined as

750

751 .. math ::

752 H_\nu^{(1)} = J_\nu(z) + iY_\nu(z),

753

754 where $J_\nu(z)$ is the Bessel function of the first kind, and

755 $Y_\nu(z)$ is the Bessel function of the second kind.

756

757 It is a solution to Bessel's equation.

758

759 Examples

760 ========

761

762 >>> from sympy import hankel1

763 >>> from sympy.abc import z, n

764 >>> hankel1(n, z).diff(z)

765 hankel1(n - 1, z)/2 - hankel1(n + 1, z)/2

766

767 See Also

768 ========

769

770 hankel2, besselj, bessely

771

772 References

773 ==========

774

775 .. [1] https://functions.wolfram.com/Bessel-TypeFunctions/HankelH1/

776

777 """

778

779 _a = S.One

780 _b = S.One

781

782 def _eval_conjugate(self):

783 z = self.argument

784 if z.is_extended_negative is False:

785 return hankel2(self.order.conjugate(), z.conjugate())

786

787

788class hankel2(BesselBase):

789 r"""

790 Hankel function of the second kind.

791

792 Explanation

793 ===========

794

795 This function is defined as

796

797 .. math ::

798 H_\nu^{(2)} = J_\nu(z) - iY_\nu(z),

799

800 where $J_\nu(z)$ is the Bessel function of the first kind, and

801 $Y_\nu(z)$ is the Bessel function of the second kind.

802

803 It is a solution to Bessel's equation, and linearly independent from

804 $H_\nu^{(1)}$.

805

806 Examples

807 ========

808

809 >>> from sympy import hankel2

810 >>> from sympy.abc import z, n

811 >>> hankel2(n, z).diff(z)

812 hankel2(n - 1, z)/2 - hankel2(n + 1, z)/2

813

814 See Also

815 ========

816

817 hankel1, besselj, bessely

818

819 References

820 ==========

821

822 .. [1] https://functions.wolfram.com/Bessel-TypeFunctions/HankelH2/

823

824 """

825

826 _a = S.One

827 _b = S.One

828

829 def _eval_conjugate(self):

830 z = self.argument

831 if z.is_extended_negative is False:

832 return hankel1(self.order.conjugate(), z.conjugate())

833

834

835def assume_integer_order(fn):

836 @wraps(fn)

837 def g(self, nu, z):

838 if nu.is_integer:

839 return fn(self, nu, z)

840 return g

841

842

843class SphericalBesselBase(BesselBase):

844 """

845 Base class for spherical Bessel functions.

846

847 These are thin wrappers around ordinary Bessel functions,

848 since spherical Bessel functions differ from the ordinary

849 ones just by a slight change in order.

850

851 To use this class, define the ``_eval_evalf()`` and ``_expand()`` methods.

852

853 """

854

855 def _expand(self, **hints):

856 """ Expand self into a polynomial. Nu is guaranteed to be Integer. """

857 raise NotImplementedError('expansion')

858

859 def _eval_expand_func(self, **hints):

860 if self.order.is_Integer:

861 return self._expand(**hints)

862 return self

863

864 def fdiff(self, argindex=2):

865 if argindex != 2:

866 raise ArgumentIndexError(self, argindex)

867 return self.__class__(self.order - 1, self.argument) - \

868 self * (self.order + 1)/self.argument

869

870

871def _jn(n, z):

872 return (spherical_bessel_fn(n, z)*sin(z) +

873 S.NegativeOne**(n + 1)*spherical_bessel_fn(-n - 1, z)*cos(z))

874

875

876def _yn(n, z):

877 # (-1)**(n + 1) * _jn(-n - 1, z)

878 return (S.NegativeOne**(n + 1) * spherical_bessel_fn(-n - 1, z)*sin(z) -

879 spherical_bessel_fn(n, z)*cos(z))

880

881

882class jn(SphericalBesselBase):

883 r"""

884 Spherical Bessel function of the first kind.

885

886 Explanation

887 ===========

888

889 This function is a solution to the spherical Bessel equation

890

891 .. math ::

892 z^2 \frac{\mathrm{d}^2 w}{\mathrm{d}z^2}

893 + 2z \frac{\mathrm{d}w}{\mathrm{d}z} + (z^2 - \nu(\nu + 1)) w = 0.

894

895 It can be defined as

896

897 .. math ::

898 j_\nu(z) = \sqrt{\frac{\pi}{2z}} J_{\nu + \frac{1}{2}}(z),

899

900 where $J_\nu(z)$ is the Bessel function of the first kind.

901

902 The spherical Bessel functions of integral order are

903 calculated using the formula:

904

905 .. math:: j_n(z) = f_n(z) \sin{z} + (-1)^{n+1} f_{-n-1}(z) \cos{z},

906

907 where the coefficients $f_n(z)$ are available as

908 :func:`sympy.polys.orthopolys.spherical_bessel_fn`.

909

910 Examples

911 ========

912

913 >>> from sympy import Symbol, jn, sin, cos, expand_func, besselj, bessely

914 >>> z = Symbol("z")

915 >>> nu = Symbol("nu", integer=True)

916 >>> print(expand_func(jn(0, z)))

917 sin(z)/z

918 >>> expand_func(jn(1, z)) == sin(z)/z**2 - cos(z)/z

919 True

920 >>> expand_func(jn(3, z))

921 (-6/z**2 + 15/z**4)*sin(z) + (1/z - 15/z**3)*cos(z)

922 >>> jn(nu, z).rewrite(besselj)

923 sqrt(2)*sqrt(pi)*sqrt(1/z)*besselj(nu + 1/2, z)/2

924 >>> jn(nu, z).rewrite(bessely)

925 (-1)**nu*sqrt(2)*sqrt(pi)*sqrt(1/z)*bessely(-nu - 1/2, z)/2

926 >>> jn(2, 5.2+0.3j).evalf(20)

927 0.099419756723640344491 - 0.054525080242173562897*I

928

929 See Also

930 ========

931

932 besselj, bessely, besselk, yn

933

934 References

935 ==========

936

937 .. [1] https://dlmf.nist.gov/10.47

938

939 """

940 @classmethod

941 def eval(cls, nu, z):

942 if z.is_zero:

943 if nu.is_zero:

944 return S.One

945 elif nu.is_integer:

946 if nu.is_positive:

947 return S.Zero

948 else:

949 return S.ComplexInfinity

950 if z in (S.NegativeInfinity, S.Infinity):

951 return S.Zero

952

953 def _eval_rewrite_as_besselj(self, nu, z, **kwargs):

954 return sqrt(pi/(2*z)) * besselj(nu + S.Half, z)

955

956 def _eval_rewrite_as_bessely(self, nu, z, **kwargs):

957 return S.NegativeOne**nu * sqrt(pi/(2*z)) * bessely(-nu - S.Half, z)

958

959 def _eval_rewrite_as_yn(self, nu, z, **kwargs):

960 return S.NegativeOne**(nu) * yn(-nu - 1, z)

961

962 def _expand(self, **hints):

963 return _jn(self.order, self.argument)

964

965 def _eval_evalf(self, prec):

966 if self.order.is_Integer:

967 return self.rewrite(besselj)._eval_evalf(prec)

968

969

970class yn(SphericalBesselBase):

971 r"""

972 Spherical Bessel function of the second kind.

973

974 Explanation

975 ===========

976

977 This function is another solution to the spherical Bessel equation, and

978 linearly independent from $j_n$. It can be defined as

979

980 .. math ::

981 y_\nu(z) = \sqrt{\frac{\pi}{2z}} Y_{\nu + \frac{1}{2}}(z),

982

983 where $Y_\nu(z)$ is the Bessel function of the second kind.

984

985 For integral orders $n$, $y_n$ is calculated using the formula:

986

987 .. math:: y_n(z) = (-1)^{n+1} j_{-n-1}(z)

988

989 Examples

990 ========

991

992 >>> from sympy import Symbol, yn, sin, cos, expand_func, besselj, bessely

993 >>> z = Symbol("z")

994 >>> nu = Symbol("nu", integer=True)

995 >>> print(expand_func(yn(0, z)))

996 -cos(z)/z

997 >>> expand_func(yn(1, z)) == -cos(z)/z**2-sin(z)/z

998 True

999 >>> yn(nu, z).rewrite(besselj)

1000 (-1)**(nu + 1)*sqrt(2)*sqrt(pi)*sqrt(1/z)*besselj(-nu - 1/2, z)/2

1001 >>> yn(nu, z).rewrite(bessely)

1002 sqrt(2)*sqrt(pi)*sqrt(1/z)*bessely(nu + 1/2, z)/2

1003 >>> yn(2, 5.2+0.3j).evalf(20)

1004 0.18525034196069722536 + 0.014895573969924817587*I

1005

1006 See Also

1007 ========

1008

1009 besselj, bessely, besselk, jn

1010

1011 References

1012 ==========

1013

1014 .. [1] https://dlmf.nist.gov/10.47

1015

1016 """

1017 @assume_integer_order

1018 def _eval_rewrite_as_besselj(self, nu, z, **kwargs):

1019 return S.NegativeOne**(nu+1) * sqrt(pi/(2*z)) * besselj(-nu - S.Half, z)

1020

1021 @assume_integer_order

1022 def _eval_rewrite_as_bessely(self, nu, z, **kwargs):

1023 return sqrt(pi/(2*z)) * bessely(nu + S.Half, z)

1024

1025 def _eval_rewrite_as_jn(self, nu, z, **kwargs):

1026 return S.NegativeOne**(nu + 1) * jn(-nu - 1, z)

1027

1028 def _expand(self, **hints):

1029 return _yn(self.order, self.argument)

1030

1031 def _eval_evalf(self, prec):

1032 if self.order.is_Integer:

1033 return self.rewrite(bessely)._eval_evalf(prec)

1034

1035

1036class SphericalHankelBase(SphericalBesselBase):

1037

1038 @assume_integer_order

1039 def _eval_rewrite_as_besselj(self, nu, z, **kwargs):

1040 # jn +- I*yn

1041 # jn as beeselj: sqrt(pi/(2*z)) * besselj(nu + S.Half, z)

1042 # yn as besselj: (-1)**(nu+1) * sqrt(pi/(2*z)) * besselj(-nu - S.Half, z)

1043 hks = self._hankel_kind_sign

1044 return sqrt(pi/(2*z))*(besselj(nu + S.Half, z) +

1045 hks*I*S.NegativeOne**(nu+1)*besselj(-nu - S.Half, z))

1046

1047 @assume_integer_order

1048 def _eval_rewrite_as_bessely(self, nu, z, **kwargs):

1049 # jn +- I*yn

1050 # jn as bessely: (-1)**nu * sqrt(pi/(2*z)) * bessely(-nu - S.Half, z)

1051 # yn as bessely: sqrt(pi/(2*z)) * bessely(nu + S.Half, z)

1052 hks = self._hankel_kind_sign

1053 return sqrt(pi/(2*z))*(S.NegativeOne**nu*bessely(-nu - S.Half, z) +

1054 hks*I*bessely(nu + S.Half, z))

1055

1056 def _eval_rewrite_as_yn(self, nu, z, **kwargs):

1057 hks = self._hankel_kind_sign

1058 return jn(nu, z).rewrite(yn) + hks*I*yn(nu, z)

1059

1060 def _eval_rewrite_as_jn(self, nu, z, **kwargs):

1061 hks = self._hankel_kind_sign

1062 return jn(nu, z) + hks*I*yn(nu, z).rewrite(jn)

1063

1064 def _eval_expand_func(self, **hints):

1065 if self.order.is_Integer:

1066 return self._expand(**hints)

1067 else:

1068 nu = self.order

1069 z = self.argument

1070 hks = self._hankel_kind_sign

1071 return jn(nu, z) + hks*I*yn(nu, z)

1072

1073 def _expand(self, **hints):

1074 n = self.order

1075 z = self.argument

1076 hks = self._hankel_kind_sign

1077

1078 # fully expanded version

1079 # return ((fn(n, z) * sin(z) +

1080 # (-1)**(n + 1) * fn(-n - 1, z) * cos(z)) + # jn

1081 # (hks * I * (-1)**(n + 1) *

1082 # (fn(-n - 1, z) * hk * I * sin(z) +

1083 # (-1)**(-n) * fn(n, z) * I * cos(z))) # +-I*yn

1084 # )

1085

1086 return (_jn(n, z) + hks*I*_yn(n, z)).expand()

1087

1088 def _eval_evalf(self, prec):

1089 if self.order.is_Integer:

1090 return self.rewrite(besselj)._eval_evalf(prec)

1091

1092

1093class hn1(SphericalHankelBase):

1094 r"""

1095 Spherical Hankel function of the first kind.

1096

1097 Explanation

1098 ===========

1099

1100 This function is defined as

1101

1102 .. math:: h_\nu^(1)(z) = j_\nu(z) + i y_\nu(z),

1103

1104 where $j_\nu(z)$ and $y_\nu(z)$ are the spherical

1105 Bessel function of the first and second kinds.

1106

1107 For integral orders $n$, $h_n^(1)$ is calculated using the formula:

1108

1109 .. math:: h_n^(1)(z) = j_{n}(z) + i (-1)^{n+1} j_{-n-1}(z)

1110

1111 Examples

1112 ========

1113

1114 >>> from sympy import Symbol, hn1, hankel1, expand_func, yn, jn

1115 >>> z = Symbol("z")

1116 >>> nu = Symbol("nu", integer=True)

1117 >>> print(expand_func(hn1(nu, z)))

1118 jn(nu, z) + I*yn(nu, z)

1119 >>> print(expand_func(hn1(0, z)))

1120 sin(z)/z - I*cos(z)/z

1121 >>> print(expand_func(hn1(1, z)))

1122 -I*sin(z)/z - cos(z)/z + sin(z)/z**2 - I*cos(z)/z**2

1123 >>> hn1(nu, z).rewrite(jn)

1124 (-1)**(nu + 1)*I*jn(-nu - 1, z) + jn(nu, z)

1125 >>> hn1(nu, z).rewrite(yn)

1126 (-1)**nu*yn(-nu - 1, z) + I*yn(nu, z)

1127 >>> hn1(nu, z).rewrite(hankel1)

1128 sqrt(2)*sqrt(pi)*sqrt(1/z)*hankel1(nu, z)/2

1129

1130 See Also

1131 ========

1132

1133 hn2, jn, yn, hankel1, hankel2

1134

1135 References

1136 ==========

1137

1138 .. [1] https://dlmf.nist.gov/10.47

1139

1140 """

1141

1142 _hankel_kind_sign = S.One

1143

1144 @assume_integer_order

1145 def _eval_rewrite_as_hankel1(self, nu, z, **kwargs):

1146 return sqrt(pi/(2*z))*hankel1(nu, z)

1147

1148

1149class hn2(SphericalHankelBase):

1150 r"""

1151 Spherical Hankel function of the second kind.

1152

1153 Explanation

1154 ===========

1155

1156 This function is defined as

1157

1158 .. math:: h_\nu^(2)(z) = j_\nu(z) - i y_\nu(z),

1159

1160 where $j_\nu(z)$ and $y_\nu(z)$ are the spherical

1161 Bessel function of the first and second kinds.

1162

1163 For integral orders $n$, $h_n^(2)$ is calculated using the formula:

1164

1165 .. math:: h_n^(2)(z) = j_{n} - i (-1)^{n+1} j_{-n-1}(z)

1166

1167 Examples

1168 ========

1169

1170 >>> from sympy import Symbol, hn2, hankel2, expand_func, jn, yn

1171 >>> z = Symbol("z")

1172 >>> nu = Symbol("nu", integer=True)

1173 >>> print(expand_func(hn2(nu, z)))

1174 jn(nu, z) - I*yn(nu, z)

1175 >>> print(expand_func(hn2(0, z)))

1176 sin(z)/z + I*cos(z)/z

1177 >>> print(expand_func(hn2(1, z)))

1178 I*sin(z)/z - cos(z)/z + sin(z)/z**2 + I*cos(z)/z**2

1179 >>> hn2(nu, z).rewrite(hankel2)

1180 sqrt(2)*sqrt(pi)*sqrt(1/z)*hankel2(nu, z)/2

1181 >>> hn2(nu, z).rewrite(jn)

1182 -(-1)**(nu + 1)*I*jn(-nu - 1, z) + jn(nu, z)

1183 >>> hn2(nu, z).rewrite(yn)

1184 (-1)**nu*yn(-nu - 1, z) - I*yn(nu, z)

1185

1186 See Also

1187 ========

1188

1189 hn1, jn, yn, hankel1, hankel2

1190

1191 References

1192 ==========

1193

1194 .. [1] https://dlmf.nist.gov/10.47

1195

1196 """

1197

1198 _hankel_kind_sign = -S.One

1199

1200 @assume_integer_order

1201 def _eval_rewrite_as_hankel2(self, nu, z, **kwargs):

1202 return sqrt(pi/(2*z))*hankel2(nu, z)

1203

1204

1205def jn_zeros(n, k, method="sympy", dps=15):

1206 """

1207 Zeros of the spherical Bessel function of the first kind.

1208

1209 Explanation

1210 ===========

1211

1212 This returns an array of zeros of $jn$ up to the $k$-th zero.

1213

1214 * method = "sympy": uses `mpmath.besseljzero

1215 <https://mpmath.org/doc/current/functions/bessel.html#mpmath.besseljzero>`_

1216 * method = "scipy": uses the

1217 `SciPy's sph_jn <https://docs.scipy.org/doc/scipy/reference/generated/scipy.special.jn_zeros.html>`_

1218 and

1219 `newton <https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.newton.html>`_

1220 to find all

1221 roots, which is faster than computing the zeros using a general

1222 numerical solver, but it requires SciPy and only works with low

1223 precision floating point numbers. (The function used with

1224 method="sympy" is a recent addition to mpmath; before that a general

1225 solver was used.)

1226

1227 Examples

1228 ========

1229

1230 >>> from sympy import jn_zeros

1231 >>> jn_zeros(2, 4, dps=5)

1232 [5.7635, 9.095, 12.323, 15.515]

1233

1234 See Also

1235 ========

1236

1237 jn, yn, besselj, besselk, bessely

1238

1239 Parameters

1240 ==========

1241

1242 n : integer

1243 order of Bessel function

1244

1245 k : integer

1246 number of zeros to return

1247

1248

1249 """

1250 from math import pi as math_pi

1251

1252 if method == "sympy":

1253 from mpmath import besseljzero

1254 from mpmath.libmp.libmpf import dps_to_prec

1255 prec = dps_to_prec(dps)

1256 return [Expr._from_mpmath(besseljzero(S(n + 0.5)._to_mpmath(prec),

1257 int(l)), prec)

1258 for l in range(1, k + 1)]

1259 elif method == "scipy":

1260 from scipy.optimize import newton

1261 try:

1262 from scipy.special import spherical_jn

1263 f = lambda x: spherical_jn(n, x)

1264 except ImportError:

1265 from scipy.special import sph_jn

1266 f = lambda x: sph_jn(n, x)[0][-1]

1267 else:

1268 raise NotImplementedError("Unknown method.")

1269

1270 def solver(f, x):

1271 if method == "scipy":

1272 root = newton(f, x)

1273 else:

1274 raise NotImplementedError("Unknown method.")

1275 return root

1276

1277 # we need to approximate the position of the first root:

1278 root = n + math_pi

1279 # determine the first root exactly:

1280 root = solver(f, root)

1281 roots = [root]

1282 for i in range(k - 1):

1283 # estimate the position of the next root using the last root + pi:

1284 root = solver(f, root + math_pi)

1285 roots.append(root)

1286 return roots

1287

1288

1289class AiryBase(Function):

1290 """

1291 Abstract base class for Airy functions.

1292

1293 This class is meant to reduce code duplication.

1294

1295 """

1296

1297 def _eval_conjugate(self):

1298 return self.func(self.args[0].conjugate())

1299

1300 def _eval_is_extended_real(self):

1301 return self.args[0].is_extended_real

1302

1303 def as_real_imag(self, deep=True, **hints):

1304 z = self.args[0]

1305 zc = z.conjugate()

1306 f = self.func

1307 u = (f(z)+f(zc))/2

1308 v = I*(f(zc)-f(z))/2

1309 return u, v

1310

1311 def _eval_expand_complex(self, deep=True, **hints):

1312 re_part, im_part = self.as_real_imag(deep=deep, **hints)

1313 return re_part + im_part*I

1314

1315

1316class airyai(AiryBase):

1317 r"""

1318 The Airy function $\operatorname{Ai}$ of the first kind.

1319

1320 Explanation

1321 ===========

1322

1323 The Airy function $\operatorname{Ai}(z)$ is defined to be the function

1324 satisfying Airy's differential equation

1325

1326 .. math::

1327 \frac{\mathrm{d}^2 w(z)}{\mathrm{d}z^2} - z w(z) = 0.

1328

1329 Equivalently, for real $z$

1330

1331 .. math::

1332 \operatorname{Ai}(z) := \frac{1}{\pi}

1333 \int_0^\infty \cos\left(\frac{t^3}{3} + z t\right) \mathrm{d}t.

1334

1335 Examples

1336 ========

1337

1338 Create an Airy function object:

1339

1340 >>> from sympy import airyai

1341 >>> from sympy.abc import z

1342

1343 >>> airyai(z)

1344 airyai(z)

1345

1346 Several special values are known:

1347

1348 >>> airyai(0)

1349 3**(1/3)/(3*gamma(2/3))

1350 >>> from sympy import oo

1351 >>> airyai(oo)

1352 0

1353 >>> airyai(-oo)

1354 0

1355

1356 The Airy function obeys the mirror symmetry:

1357

1358 >>> from sympy import conjugate

1359 >>> conjugate(airyai(z))

1360 airyai(conjugate(z))

1361

1362 Differentiation with respect to $z$ is supported:

1363

1364 >>> from sympy import diff

1365 >>> diff(airyai(z), z)

1366 airyaiprime(z)

1367 >>> diff(airyai(z), z, 2)

1368 z*airyai(z)

1369

1370 Series expansion is also supported:

1371

1372 >>> from sympy import series

1373 >>> series(airyai(z), z, 0, 3)

1374 3**(5/6)*gamma(1/3)/(6*pi) - 3**(1/6)*z*gamma(2/3)/(2*pi) + O(z**3)

1375

1376 We can numerically evaluate the Airy function to arbitrary precision

1377 on the whole complex plane:

1378

1379 >>> airyai(-2).evalf(50)

1380 0.22740742820168557599192443603787379946077222541710

1381

1382 Rewrite $\operatorname{Ai}(z)$ in terms of hypergeometric functions:

1383

1384 >>> from sympy import hyper

1385 >>> airyai(z).rewrite(hyper)

1386 -3**(2/3)*z*hyper((), (4/3,), z**3/9)/(3*gamma(1/3)) + 3**(1/3)*hyper((), (2/3,), z**3/9)/(3*gamma(2/3))

1387

1388 See Also

1389 ========

1390

1391 airybi: Airy function of the second kind.

1392 airyaiprime: Derivative of the Airy function of the first kind.

1393 airybiprime: Derivative of the Airy function of the second kind.

1394

1395 References

1396 ==========

1397

1398 .. [1] https://en.wikipedia.org/wiki/Airy_function

1399 .. [2] https://dlmf.nist.gov/9

1400 .. [3] https://encyclopediaofmath.org/wiki/Airy_functions

1401 .. [4] https://mathworld.wolfram.com/AiryFunctions.html

1402

1403 """

1404

1405 nargs = 1

1406 unbranched = True

1407

1408 @classmethod

1409 def eval(cls, arg):

1410 if arg.is_Number:

1411 if arg is S.NaN:

1412 return S.NaN

1413 elif arg is S.Infinity:

1414 return S.Zero

1415 elif arg is S.NegativeInfinity:

1416 return S.Zero

1417 elif arg.is_zero:

1418 return S.One / (3**Rational(2, 3) * gamma(Rational(2, 3)))

1419 if arg.is_zero:

1420 return S.One / (3**Rational(2, 3) * gamma(Rational(2, 3)))

1421

1422 def fdiff(self, argindex=1):

1423 if argindex == 1:

1424 return airyaiprime(self.args[0])

1425 else:

1426 raise ArgumentIndexError(self, argindex)

1427

1428 @staticmethod

1429 @cacheit

1430 def taylor_term(n, x, *previous_terms):

1431 if n < 0:

1432 return S.Zero

1433 else:

1434 x = sympify(x)

1435 if len(previous_terms) > 1:

1436 p = previous_terms[-1]

1437 return ((cbrt(3)*x)**(-n)*(cbrt(3)*x)**(n + 1)*sin(pi*(n*Rational(2, 3) + Rational(4, 3)))*factorial(n) *

1438 gamma(n/3 + Rational(2, 3))/(sin(pi*(n*Rational(2, 3) + Rational(2, 3)))*factorial(n + 1)*gamma(n/3 + Rational(1, 3))) * p)

1439 else:

1440 return (S.One/(3**Rational(2, 3)*pi) * gamma((n+S.One)/S(3)) * sin(Rational(2, 3)*pi*(n+S.One)) /

1441 factorial(n) * (cbrt(3)*x)**n)

1442

1443 def _eval_rewrite_as_besselj(self, z, **kwargs):

1444 ot = Rational(1, 3)

1445 tt = Rational(2, 3)

1446 a = Pow(-z, Rational(3, 2))

1447 if re(z).is_negative:

1448 return ot*sqrt(-z) * (besselj(-ot, tt*a) + besselj(ot, tt*a))

1449

1450 def _eval_rewrite_as_besseli(self, z, **kwargs):

1451 ot = Rational(1, 3)

1452 tt = Rational(2, 3)

1453 a = Pow(z, Rational(3, 2))

1454 if re(z).is_positive:

1455 return ot*sqrt(z) * (besseli(-ot, tt*a) - besseli(ot, tt*a))

1456 else:

1457 return ot*(Pow(a, ot)*besseli(-ot, tt*a) - z*Pow(a, -ot)*besseli(ot, tt*a))

1458

1459 def _eval_rewrite_as_hyper(self, z, **kwargs):

1460 pf1 = S.One / (3**Rational(2, 3)*gamma(Rational(2, 3)))

1461 pf2 = z / (root(3, 3)*gamma(Rational(1, 3)))

1462 return pf1 * hyper([], [Rational(2, 3)], z**3/9) - pf2 * hyper([], [Rational(4, 3)], z**3/9)

1463

1464 def _eval_expand_func(self, **hints):

1465 arg = self.args[0]

1466 symbs = arg.free_symbols

1467

1468 if len(symbs) == 1:

1469 z = symbs.pop()

1470 c = Wild("c", exclude=[z])

1471 d = Wild("d", exclude=[z])

1472 m = Wild("m", exclude=[z])

1473 n = Wild("n", exclude=[z])

1474 M = arg.match(c*(d*z**n)**m)

1475 if M is not None:

1476 m = M[m]

1477 # The transformation is given by 03.05.16.0001.01

1478 # https://functions.wolfram.com/Bessel-TypeFunctions/AiryAi/16/01/01/0001/

1479 if (3*m).is_integer:

1480 c = M[c]

1481 d = M[d]

1482 n = M[n]

1483 pf = (d * z**n)**m / (d**m * z**(m*n))

1484 newarg = c * d**m * z**(m*n)

1485 return S.Half * ((pf + S.One)*airyai(newarg) - (pf - S.One)/sqrt(3)*airybi(newarg))

1486

1487

1488class airybi(AiryBase):

1489 r"""

1490 The Airy function $\operatorname{Bi}$ of the second kind.

1491

1492 Explanation

1493 ===========

1494

1495 The Airy function $\operatorname{Bi}(z)$ is defined to be the function

1496 satisfying Airy's differential equation

1497

1498 .. math::

1499 \frac{\mathrm{d}^2 w(z)}{\mathrm{d}z^2} - z w(z) = 0.

1500

1501 Equivalently, for real $z$

1502

1503 .. math::

1504 \operatorname{Bi}(z) := \frac{1}{\pi}

1505 \int_0^\infty

1506 \exp\left(-\frac{t^3}{3} + z t\right)

1507 + \sin\left(\frac{t^3}{3} + z t\right) \mathrm{d}t.

1508

1509 Examples

1510 ========

1511

1512 Create an Airy function object:

1513

1514 >>> from sympy import airybi

1515 >>> from sympy.abc import z

1516

1517 >>> airybi(z)

1518 airybi(z)

1519

1520 Several special values are known:

1521

1522 >>> airybi(0)

1523 3**(5/6)/(3*gamma(2/3))

1524 >>> from sympy import oo

1525 >>> airybi(oo)

1526 oo

1527 >>> airybi(-oo)

1528 0

1529

1530 The Airy function obeys the mirror symmetry:

1531

1532 >>> from sympy import conjugate

1533 >>> conjugate(airybi(z))

1534 airybi(conjugate(z))

1535

1536 Differentiation with respect to $z$ is supported:

1537

1538 >>> from sympy import diff

1539 >>> diff(airybi(z), z)

1540 airybiprime(z)

1541 >>> diff(airybi(z), z, 2)

1542 z*airybi(z)

1543

1544 Series expansion is also supported:

1545

1546 >>> from sympy import series

1547 >>> series(airybi(z), z, 0, 3)

1548 3**(1/3)*gamma(1/3)/(2*pi) + 3**(2/3)*z*gamma(2/3)/(2*pi) + O(z**3)

1549

1550 We can numerically evaluate the Airy function to arbitrary precision

1551 on the whole complex plane:

1552

1553 >>> airybi(-2).evalf(50)

1554 -0.41230258795639848808323405461146104203453483447240

1555

1556 Rewrite $\operatorname{Bi}(z)$ in terms of hypergeometric functions:

1557

1558 >>> from sympy import hyper

1559 >>> airybi(z).rewrite(hyper)

1560 3**(1/6)*z*hyper((), (4/3,), z**3/9)/gamma(1/3) + 3**(5/6)*hyper((), (2/3,), z**3/9)/(3*gamma(2/3))

1561

1562 See Also

1563 ========

1564

1565 airyai: Airy function of the first kind.

1566 airyaiprime: Derivative of the Airy function of the first kind.

1567 airybiprime: Derivative of the Airy function of the second kind.

1568

1569 References

1570 ==========

1571

1572 .. [1] https://en.wikipedia.org/wiki/Airy_function

1573 .. [2] https://dlmf.nist.gov/9

1574 .. [3] https://encyclopediaofmath.org/wiki/Airy_functions

1575 .. [4] https://mathworld.wolfram.com/AiryFunctions.html

1576

1577 """

1578

1579 nargs = 1

1580 unbranched = True

1581

1582 @classmethod

1583 def eval(cls, arg):

1584 if arg.is_Number:

1585 if arg is S.NaN:

1586 return S.NaN

1587 elif arg is S.Infinity:

1588 return S.Infinity

1589 elif arg is S.NegativeInfinity:

1590 return S.Zero

1591 elif arg.is_zero:

1592 return S.One / (3**Rational(1, 6) * gamma(Rational(2, 3)))

1593

1594 if arg.is_zero:

1595 return S.One / (3**Rational(1, 6) * gamma(Rational(2, 3)))

1596

1597 def fdiff(self, argindex=1):

1598 if argindex == 1:

1599 return airybiprime(self.args[0])

1600 else:

1601 raise ArgumentIndexError(self, argindex)

1602

1603 @staticmethod

1604 @cacheit

1605 def taylor_term(n, x, *previous_terms):

1606 if n < 0:

1607 return S.Zero

1608 else:

1609 x = sympify(x)

1610 if len(previous_terms) > 1:

1611 p = previous_terms[-1]

1612 return (cbrt(3)*x * Abs(sin(Rational(2, 3)*pi*(n + S.One))) * factorial((n - S.One)/S(3)) /

1613 ((n + S.One) * Abs(cos(Rational(2, 3)*pi*(n + S.Half))) * factorial((n - 2)/S(3))) * p)

1614 else:

1615 return (S.One/(root(3, 6)*pi) * gamma((n + S.One)/S(3)) * Abs(sin(Rational(2, 3)*pi*(n + S.One))) /

1616 factorial(n) * (cbrt(3)*x)**n)

1617

1618 def _eval_rewrite_as_besselj(self, z, **kwargs):

1619 ot = Rational(1, 3)

1620 tt = Rational(2, 3)

1621 a = Pow(-z, Rational(3, 2))

1622 if re(z).is_negative:

1623 return sqrt(-z/3) * (besselj(-ot, tt*a) - besselj(ot, tt*a))

1624

1625 def _eval_rewrite_as_besseli(self, z, **kwargs):

1626 ot = Rational(1, 3)

1627 tt = Rational(2, 3)

1628 a = Pow(z, Rational(3, 2))

1629 if re(z).is_positive:

1630 return sqrt(z)/sqrt(3) * (besseli(-ot, tt*a) + besseli(ot, tt*a))

1631 else:

1632 b = Pow(a, ot)

1633 c = Pow(a, -ot)

1634 return sqrt(ot)*(b*besseli(-ot, tt*a) + z*c*besseli(ot, tt*a))

1635

1636 def _eval_rewrite_as_hyper(self, z, **kwargs):

1637 pf1 = S.One / (root(3, 6)*gamma(Rational(2, 3)))

1638 pf2 = z*root(3, 6) / gamma(Rational(1, 3))

1639 return pf1 * hyper([], [Rational(2, 3)], z**3/9) + pf2 * hyper([], [Rational(4, 3)], z**3/9)

1640

1641 def _eval_expand_func(self, **hints):

1642 arg = self.args[0]

1643 symbs = arg.free_symbols

1644

1645 if len(symbs) == 1:

1646 z = symbs.pop()

1647 c = Wild("c", exclude=[z])

1648 d = Wild("d", exclude=[z])

1649 m = Wild("m", exclude=[z])

1650 n = Wild("n", exclude=[z])

1651 M = arg.match(c*(d*z**n)**m)

1652 if M is not None:

1653 m = M[m]

1654 # The transformation is given by 03.06.16.0001.01

1655 # https://functions.wolfram.com/Bessel-TypeFunctions/AiryBi/16/01/01/0001/

1656 if (3*m).is_integer:

1657 c = M[c]

1658 d = M[d]

1659 n = M[n]

1660 pf = (d * z**n)**m / (d**m * z**(m*n))

1661 newarg = c * d**m * z**(m*n)

1662 return S.Half * (sqrt(3)*(S.One - pf)*airyai(newarg) + (S.One + pf)*airybi(newarg))

1663

1664

1665class airyaiprime(AiryBase):

1666 r"""

1667 The derivative $\operatorname{Ai}^\prime$ of the Airy function of the first

1668 kind.

1669

1670 Explanation

1671 ===========

1672

1673 The Airy function $\operatorname{Ai}^\prime(z)$ is defined to be the

1674 function

1675

1676 .. math::

1677 \operatorname{Ai}^\prime(z) := \frac{\mathrm{d} \operatorname{Ai}(z)}{\mathrm{d} z}.

1678

1679 Examples

1680 ========

1681

1682 Create an Airy function object:

1683

1684 >>> from sympy import airyaiprime

1685 >>> from sympy.abc import z

1686

1687 >>> airyaiprime(z)

1688 airyaiprime(z)

1689

1690 Several special values are known:

1691

1692 >>> airyaiprime(0)

1693 -3**(2/3)/(3*gamma(1/3))

1694 >>> from sympy import oo

1695 >>> airyaiprime(oo)

1696 0

1697

1698 The Airy function obeys the mirror symmetry:

1699

1700 >>> from sympy import conjugate

1701 >>> conjugate(airyaiprime(z))

1702 airyaiprime(conjugate(z))

1703

1704 Differentiation with respect to $z$ is supported:

1705

1706 >>> from sympy import diff

1707 >>> diff(airyaiprime(z), z)

1708 z*airyai(z)

1709 >>> diff(airyaiprime(z), z, 2)

1710 z*airyaiprime(z) + airyai(z)

1711

1712 Series expansion is also supported:

1713

1714 >>> from sympy import series

1715 >>> series(airyaiprime(z), z, 0, 3)

1716 -3**(2/3)/(3*gamma(1/3)) + 3**(1/3)*z**2/(6*gamma(2/3)) + O(z**3)

1717

1718 We can numerically evaluate the Airy function to arbitrary precision

1719 on the whole complex plane:

1720

1721 >>> airyaiprime(-2).evalf(50)

1722 0.61825902074169104140626429133247528291577794512415

1723

1724 Rewrite $\operatorname{Ai}^\prime(z)$ in terms of hypergeometric functions:

1725

1726 >>> from sympy import hyper

1727 >>> airyaiprime(z).rewrite(hyper)

1728 3**(1/3)*z**2*hyper((), (5/3,), z**3/9)/(6*gamma(2/3)) - 3**(2/3)*hyper((), (1/3,), z**3/9)/(3*gamma(1/3))

1729

1730 See Also

1731 ========

1732

1733 airyai: Airy function of the first kind.

1734 airybi: Airy function of the second kind.

1735 airybiprime: Derivative of the Airy function of the second kind.

1736

1737 References

1738 ==========

1739

1740 .. [1] https://en.wikipedia.org/wiki/Airy_function

1741 .. [2] https://dlmf.nist.gov/9

1742 .. [3] https://encyclopediaofmath.org/wiki/Airy_functions

1743 .. [4] https://mathworld.wolfram.com/AiryFunctions.html

1744

1745 """

1746

1747 nargs = 1

1748 unbranched = True

1749

1750 @classmethod

1751 def eval(cls, arg):

1752 if arg.is_Number:

1753 if arg is S.NaN:

1754 return S.NaN

1755 elif arg is S.Infinity:

1756 return S.Zero

1757

1758 if arg.is_zero:

1759 return S.NegativeOne / (3**Rational(1, 3) * gamma(Rational(1, 3)))

1760

1761 def fdiff(self, argindex=1):

1762 if argindex == 1:

1763 return self.args[0]*airyai(self.args[0])

1764 else:

1765 raise ArgumentIndexError(self, argindex)

1766

1767 def _eval_evalf(self, prec):

1768 z = self.args[0]._to_mpmath(prec)

1769 with workprec(prec):

1770 res = mp.airyai(z, derivative=1)

1771 return Expr._from_mpmath(res, prec)

1772

1773 def _eval_rewrite_as_besselj(self, z, **kwargs):

1774 tt = Rational(2, 3)

1775 a = Pow(-z, Rational(3, 2))

1776 if re(z).is_negative:

1777 return z/3 * (besselj(-tt, tt*a) - besselj(tt, tt*a))

1778

1779 def _eval_rewrite_as_besseli(self, z, **kwargs):

1780 ot = Rational(1, 3)

1781 tt = Rational(2, 3)

1782 a = tt * Pow(z, Rational(3, 2))

1783 if re(z).is_positive:

1784 return z/3 * (besseli(tt, a) - besseli(-tt, a))

1785 else:

1786 a = Pow(z, Rational(3, 2))

1787 b = Pow(a, tt)

1788 c = Pow(a, -tt)

1789 return ot * (z**2*c*besseli(tt, tt*a) - b*besseli(-ot, tt*a))

1790

1791 def _eval_rewrite_as_hyper(self, z, **kwargs):

1792 pf1 = z**2 / (2*3**Rational(2, 3)*gamma(Rational(2, 3)))

1793 pf2 = 1 / (root(3, 3)*gamma(Rational(1, 3)))

1794 return pf1 * hyper([], [Rational(5, 3)], z**3/9) - pf2 * hyper([], [Rational(1, 3)], z**3/9)

1795

1796 def _eval_expand_func(self, **hints):

1797 arg = self.args[0]

1798 symbs = arg.free_symbols

1799

1800 if len(symbs) == 1:

1801 z = symbs.pop()

1802 c = Wild("c", exclude=[z])

1803 d = Wild("d", exclude=[z])

1804 m = Wild("m", exclude=[z])

1805 n = Wild("n", exclude=[z])

1806 M = arg.match(c*(d*z**n)**m)

1807 if M is not None:

1808 m = M[m]

1809 # The transformation is in principle

1810 # given by 03.07.16.0001.01 but note

1811 # that there is an error in this formula.

1812 # https://functions.wolfram.com/Bessel-TypeFunctions/AiryAiPrime/16/01/01/0001/

1813 if (3*m).is_integer:

1814 c = M[c]

1815 d = M[d]

1816 n = M[n]

1817 pf = (d**m * z**(n*m)) / (d * z**n)**m

1818 newarg = c * d**m * z**(n*m)

1819 return S.Half * ((pf + S.One)*airyaiprime(newarg) + (pf - S.One)/sqrt(3)*airybiprime(newarg))

1820

1821

1822class airybiprime(AiryBase):

1823 r"""

1824 The derivative $\operatorname{Bi}^\prime$ of the Airy function of the first

1825 kind.

1826

1827 Explanation

1828 ===========

1829

1830 The Airy function $\operatorname{Bi}^\prime(z)$ is defined to be the

1831 function

1832

1833 .. math::

1834 \operatorname{Bi}^\prime(z) := \frac{\mathrm{d} \operatorname{Bi}(z)}{\mathrm{d} z}.

1835

1836 Examples

1837 ========

1838

1839 Create an Airy function object:

1840

1841 >>> from sympy import airybiprime

1842 >>> from sympy.abc import z

1843

1844 >>> airybiprime(z)

1845 airybiprime(z)

1846

1847 Several special values are known:

1848

1849 >>> airybiprime(0)

1850 3**(1/6)/gamma(1/3)

1851 >>> from sympy import oo

1852 >>> airybiprime(oo)

1853 oo

1854 >>> airybiprime(-oo)

1855 0

1856

1857 The Airy function obeys the mirror symmetry:

1858

1859 >>> from sympy import conjugate

1860 >>> conjugate(airybiprime(z))

1861 airybiprime(conjugate(z))

1862

1863 Differentiation with respect to $z$ is supported:

1864

1865 >>> from sympy import diff

1866 >>> diff(airybiprime(z), z)

1867 z*airybi(z)

1868 >>> diff(airybiprime(z), z, 2)

1869 z*airybiprime(z) + airybi(z)

1870

1871 Series expansion is also supported:

1872

1873 >>> from sympy import series

1874 >>> series(airybiprime(z), z, 0, 3)

1875 3**(1/6)/gamma(1/3) + 3**(5/6)*z**2/(6*gamma(2/3)) + O(z**3)

1876

1877 We can numerically evaluate the Airy function to arbitrary precision

1878 on the whole complex plane:

1879

1880 >>> airybiprime(-2).evalf(50)

1881 0.27879516692116952268509756941098324140300059345163

1882

1883 Rewrite $\operatorname{Bi}^\prime(z)$ in terms of hypergeometric functions:

1884

1885 >>> from sympy import hyper

1886 >>> airybiprime(z).rewrite(hyper)

1887 3**(5/6)*z**2*hyper((), (5/3,), z**3/9)/(6*gamma(2/3)) + 3**(1/6)*hyper((), (1/3,), z**3/9)/gamma(1/3)

1888

1889 See Also

1890 ========

1891

1892 airyai: Airy function of the first kind.

1893 airybi: Airy function of the second kind.

1894 airyaiprime: Derivative of the Airy function of the first kind.

1895

1896 References

1897 ==========

1898

1899 .. [1] https://en.wikipedia.org/wiki/Airy_function

1900 .. [2] https://dlmf.nist.gov/9

1901 .. [3] https://encyclopediaofmath.org/wiki/Airy_functions

1902 .. [4] https://mathworld.wolfram.com/AiryFunctions.html

1903

1904 """

1905

1906 nargs = 1

1907 unbranched = True

1908

1909 @classmethod

1910 def eval(cls, arg):

1911 if arg.is_Number:

1912 if arg is S.NaN:

1913 return S.NaN

1914 elif arg is S.Infinity:

1915 return S.Infinity

1916 elif arg is S.NegativeInfinity:

1917 return S.Zero

1918 elif arg.is_zero:

1919 return 3**Rational(1, 6) / gamma(Rational(1, 3))

1920

1921 if arg.is_zero:

1922 return 3**Rational(1, 6) / gamma(Rational(1, 3))

1923

1924

1925 def fdiff(self, argindex=1):

1926 if argindex == 1:

1927 return self.args[0]*airybi(self.args[0])

1928 else:

1929 raise ArgumentIndexError(self, argindex)

1930

1931 def _eval_evalf(self, prec):

1932 z = self.args[0]._to_mpmath(prec)

1933 with workprec(prec):

1934 res = mp.airybi(z, derivative=1)

1935 return Expr._from_mpmath(res, prec)

1936

1937 def _eval_rewrite_as_besselj(self, z, **kwargs):

1938 tt = Rational(2, 3)

1939 a = tt * Pow(-z, Rational(3, 2))

1940 if re(z).is_negative:

1941 return -z/sqrt(3) * (besselj(-tt, a) + besselj(tt, a))

1942

1943 def _eval_rewrite_as_besseli(self, z, **kwargs):

1944 ot = Rational(1, 3)

1945 tt = Rational(2, 3)

1946 a = tt * Pow(z, Rational(3, 2))

1947 if re(z).is_positive:

1948 return z/sqrt(3) * (besseli(-tt, a) + besseli(tt, a))

1949 else:

1950 a = Pow(z, Rational(3, 2))

1951 b = Pow(a, tt)

1952 c = Pow(a, -tt)

1953 return sqrt(ot) * (b*besseli(-tt, tt*a) + z**2*c*besseli(tt, tt*a))

1954

1955 def _eval_rewrite_as_hyper(self, z, **kwargs):

1956 pf1 = z**2 / (2*root(3, 6)*gamma(Rational(2, 3)))

1957 pf2 = root(3, 6) / gamma(Rational(1, 3))

1958 return pf1 * hyper([], [Rational(5, 3)], z**3/9) + pf2 * hyper([], [Rational(1, 3)], z**3/9)

1959

1960 def _eval_expand_func(self, **hints):

1961 arg = self.args[0]

1962 symbs = arg.free_symbols

1963

1964 if len(symbs) == 1:

1965 z = symbs.pop()

1966 c = Wild("c", exclude=[z])

1967 d = Wild("d", exclude=[z])

1968 m = Wild("m", exclude=[z])

1969 n = Wild("n", exclude=[z])

1970 M = arg.match(c*(d*z**n)**m)

1971 if M is not None:

1972 m = M[m]

1973 # The transformation is in principle

1974 # given by 03.08.16.0001.01 but note

1975 # that there is an error in this formula.

1976 # https://functions.wolfram.com/Bessel-TypeFunctions/AiryBiPrime/16/01/01/0001/

1977 if (3*m).is_integer:

1978 c = M[c]

1979 d = M[d]

1980 n = M[n]

1981 pf = (d**m * z**(n*m)) / (d * z**n)**m

1982 newarg = c * d**m * z**(n*m)

1983 return S.Half * (sqrt(3)*(pf - S.One)*airyaiprime(newarg) + (pf + S.One)*airybiprime(newarg))

1984

1985

1986class marcumq(Function):

1987 r"""

1988 The Marcum Q-function.

1989

1990 Explanation

1991 ===========

1992

1993 The Marcum Q-function is defined by the meromorphic continuation of

1994

1995 .. math::

1996 Q_m(a, b) = a^{- m + 1} \int_{b}^{\infty} x^{m} e^{- \frac{a^{2}}{2} - \frac{x^{2}}{2}} I_{m - 1}\left(a x\right)\, dx

1997

1998 Examples

1999 ========

2000

2001 >>> from sympy import marcumq

2002 >>> from sympy.abc import m, a, b

2003 >>> marcumq(m, a, b)

2004 marcumq(m, a, b)

2005

2006 Special values:

2007

2008 >>> marcumq(m, 0, b)

2009 uppergamma(m, b**2/2)/gamma(m)

2010 >>> marcumq(0, 0, 0)

2011 0

2012 >>> marcumq(0, a, 0)

2013 1 - exp(-a**2/2)

2014 >>> marcumq(1, a, a)

2015 1/2 + exp(-a**2)*besseli(0, a**2)/2

2016 >>> marcumq(2, a, a)

2017 1/2 + exp(-a**2)*besseli(0, a**2)/2 + exp(-a**2)*besseli(1, a**2)

2018

2019 Differentiation with respect to $a$ and $b$ is supported:

2020

2021 >>> from sympy import diff

2022 >>> diff(marcumq(m, a, b), a)

2023 a*(-marcumq(m, a, b) + marcumq(m + 1, a, b))

2024 >>> diff(marcumq(m, a, b), b)

2025 -a**(1 - m)*b**m*exp(-a**2/2 - b**2/2)*besseli(m - 1, a*b)

2026

2027 References

2028 ==========

2029

2030 .. [1] https://en.wikipedia.org/wiki/Marcum_Q-function

2031 .. [2] https://mathworld.wolfram.com/MarcumQ-Function.html

2032

2033 """

2034

2035 @classmethod

2036 def eval(cls, m, a, b):

2037 if a is S.Zero:

2038 if m is S.Zero and b is S.Zero:

2039 return S.Zero

2040 return uppergamma(m, b**2 * S.Half) / gamma(m)

2041

2042 if m is S.Zero and b is S.Zero:

2043 return 1 - 1 / exp(a**2 * S.Half)

2044

2045 if a == b:

2046 if m is S.One:

2047 return (1 + exp(-a**2) * besseli(0, a**2))*S.Half

2048 if m == 2:

2049 return S.Half + S.Half * exp(-a**2) * besseli(0, a**2) + exp(-a**2) * besseli(1, a**2)

2050

2051 if a.is_zero:

2052 if m.is_zero and b.is_zero:

2053 return S.Zero

2054 return uppergamma(m, b**2*S.Half) / gamma(m)

2055

2056 if m.is_zero and b.is_zero:

2057 return 1 - 1 / exp(a**2*S.Half)

2058

2059 def fdiff(self, argindex=2):

2060 m, a, b = self.args

2061 if argindex == 2:

2062 return a * (-marcumq(m, a, b) + marcumq(1+m, a, b))

2063 elif argindex == 3:

2064 return (-b**m / a**(m-1)) * exp(-(a**2 + b**2)/2) * besseli(m-1, a*b)

2065 else:

2066 raise ArgumentIndexError(self, argindex)

2067

2068 def _eval_rewrite_as_Integral(self, m, a, b, **kwargs):

2069 from sympy.integrals.integrals import Integral

2070 x = kwargs.get('x', Dummy('x'))

2071 return a ** (1 - m) * \

2072 Integral(x**m * exp(-(x**2 + a**2)/2) * besseli(m-1, a*x), [x, b, S.Infinity])

2073

2074 def _eval_rewrite_as_Sum(self, m, a, b, **kwargs):

2075 from sympy.concrete.summations import Sum

2076 k = kwargs.get('k', Dummy('k'))

2077 return exp(-(a**2 + b**2) / 2) * Sum((a/b)**k * besseli(k, a*b), [k, 1-m, S.Infinity])

2078

2079 def _eval_rewrite_as_besseli(self, m, a, b, **kwargs):

2080 if a == b:

2081 if m == 1:

2082 return (1 + exp(-a**2) * besseli(0, a**2)) / 2

2083 if m.is_Integer and m >= 2:

2084 s = sum([besseli(i, a**2) for i in range(1, m)])

2085 return S.Half + exp(-a**2) * besseli(0, a**2) / 2 + exp(-a**2) * s

2086

2087 def _eval_is_zero(self):

2088 if all(arg.is_zero for arg in self.args):

2089 return True