Special Functions in Mathematica Part I: general overview Oleksandr - - PDF document

Special Functions in Mathematica

Part I: general overview

Oleksandr Pavlyk Wolfram Research Inc. pavlyk@wolfram.com

Table of Contents

Mathematica rules… Special functions… Functions in Mathematica Defining your own functions

2 Talk - Part I.nb

Mathematica rules

» Capitalized function names, e.g. Plot, Solve, Integrate » Square brackets around the arguments, i.e. Sinx, not Sin x » Curly brackets for lists and ranges » + to evaluate » Take advantage of pallettes Examples: Stand-alone commands

In[1]:= Fibonacci36 Out[1]= 14 930 352 In[2]:= Sin36 Degree Out[2]=

5 8 5 8

In[3]:= Sumx ^ k k ^ 2, k, 1, Infinity Out[3]= PolyLog2, x In[4]:= SeriesSinTanx TanSinx, x, 0, 12 Out[4]= x7

30 29 x9 756 1913 x11 75 600 Ox13

In[5]:= DSolvey ''x x yx x , y, x Out[5]= y Functionx, Π AiryAiPrimex AiryBix Π AiryAix

AiryBiPrimex AiryAix C1 AiryBix C2

Examples: Command interaction

In[6]:= sol Solve4 x 1 ^ 2 3 y ^ 2 28, 2 x 3 y 4, x, y Out[6]= x 1, y 2, x 7

2 , y 1

Talk - Part I.nb 3

In[7]:= ContourPlot4 x 1 ^ 2 3 y ^ 2 28, 2 x 3 y 4,

x, 2, 4, y, 4, 4, Epilog PointSizeLarge, Red, Pointx, y . sol

Out[7]=

2 1 1 2 3 4 4 2 2 4

4 Talk - Part I.nb

Examples: Interactive exploration

In[8]:= Manipulate

PlotExpx, Sumx ^ k k , k, 0, n, x, 3, 3, n, 0, 6, 1, SetterBar

Out[8]=

n 1 2 3 4 5 6

3 2 1 1 2 3 5 10 15 20

Talk - Part I.nb 5

Special functions…

… are particular mathematical functions which have more or less established names and notations due to their importance in mathematical analysis, functional analysis, physics, or other applications. There is no general formal definition, but the list of mathematical functions contains functions which are commonly accepted as special. In particular, elementary functions are also considered as special functions.

Comment

The above is a quote from wikipedia [link]

6 Talk - Part I.nb

… and why study them?

» to solve the problem they appeared in (physics, chemistry, statistics,

engineering, etc.)

» because they provide a language to solve other problems with

The integral

sinx x

x is not expressible in elementary functions, see Liouville’s theorem But expand the language and the integral becomes doable:

IntegrateSinx x, x SinIntegralx

Sine integral provides a language to express answers to other problems as well:

Sum 1 ^ k 2 k 1 x2 k1 2 k 1 , k, 0, SinIntegralx Integrate Expx t 1 t ^ 2 , t, 0, Infinity, Assumptions x 0 CosIntegralx Sinx 1 2 Cosx Π 2 SinIntegralx

Talk - Part I.nb 7

… and why study them?

» just for the fun of it

Example 1

Comparison of quantum and classical probability distributions for a harmonic

• scillator:

Blockn 20, Plot x2 HermiteHn, x2 2n n Π , Piecewise

• 1

2 n 1 x2 Π , x2 2 n 1, x, 3 n , 3 n

5 5 0.05 0.10 0.15 0.20

Example 2

Given a word with character counts n1, …, nk the number of anagrams with no fixed letter: DerangementsCountnvec_List : IntegrateExpx Product1n LaguerreLn, x, n, nvec, x, 0,

For “LHCPhenoNet”:

LHCPhenoNet 1211 221 1

8 Talk - Part I.nb

DerangementsCount1, 2, 1, 1, 2, 2, 1, 1 989 660 totalPermutationCount Multinomial1, 2, 1, 1, 2, 2, 1, 1 4 989 600

Talk - Part I.nb 9

Special functions in Mathematica

The chart of 180+ functions

Elliptic Functions

Hypergeometric functions

Zeta, PolyLogs Mathieu, Spheroidal Elliptic intergals Number Theory Functions , B, Ψ Statistical Functions Polynomials Holonomic q functions

Zeros : Ρk , jΝ,k Landscape of special function groups, supported in Mathematica »

10 Talk - Part I.nb

Hypergeometric function

Geometric series

• k0
• xk

1 1 x Generalized hypergeometric series

• k0
• ck xk

here ck is such that ck1 ck

1 k1 ka1 ka2kap kb1 kb2kbq .

The non-terminating series converges absolutely for p q 1, or for p q 1 and x 1. The series terminates if one ai is a non-positive integer. A function defined by an analytic continuation of the above series is denoted as

pFqa1, a2, …, ap; b1, b2, …, bq; z and in Mathematica as

HypergeometricPFQa1, a2, …, ap, b1, b2, …, bq, x

Comment

Many elementary and special functions are hypergeometric functions. Hypergeometric functions enjoy many nice properties warranting closer look in part II.

Talk - Part I.nb 11

Special functions in Mathematica

extent of the support

» Numeric evaluation » Exact evaluation » Differentiation » Series expansion » Simplification rules, functional expansion » Support in solvers (input & output)

12 Talk - Part I.nb

Special functions in Mathematica

Exact evaluation

Sequences: BernoulliB60 1 215 233 140 483 755 572 040 304 994 079 820 246 041 491 56 786 730 Binomial52, 6 20 358 520 LegendreP4, x 1 8 3 30 x2 35 x4 Evaluation at special points: ErfInfinity 1 Zeta '1 1 12 LogGlaisher Hypergeometric2F1n, n, 2 n 1, 1 22 n Gamma 1

2 n

Π Gamma1 n HermiteHn, 0 2n Π Gamma 1n

2

JacobiSNEllipticKm 2, m 1 1 1 m

Talk - Part I.nb 13

Evaluation to simpler functions: Hypergeometric2F11, 1, 4, x

• 1

2 x3 3 2 x 3 x2 2 Log1 x 4 x Log1 x 2 x2 Log1 x LegendreQ3, 2, x 1 x2 8 25 x2 15 x4 1 x22

• 15 x 1 x2 1

2 Log1 x 1 2 Log1 x Parity transformations, argument reduction make expression canonical: Erfx Erfx Sin127 Degree Cos37 ° JacobiSNu I EllipticK1 m 2 EllipticKm, m JacobiNSu, m m

14 Talk - Part I.nb

Special functions in Mathematica

Numeric evaluation: maximally extended domain

Special functions are functions complex variable. Mathematical consistency — choices of branch-cuts, and continuity at them, are consistent with functional relations between functions, and evaluation semantics. Implications of argument reduction

z2 z, because in Mathematica argument is automatically reduced to Π argz Π

Π

• Π

2 Π 2

Π Φ Π

• Π

2 Π 2

Π arg f Φ

ArgExp Φ Arg Exp Φ2 Explanation

Likewise

1 z

• 1

z , log 1 z logz, etc.

• 1

z , 1 z , Log 1 z , Logz . z 1 , , Π, Π

Talk - Part I.nb 15

ExpIntegralE1, x FunctionExpand ExpIntegralEix 1 2 Log 1 x Logx Logx Refine, x 0 ExpIntegralEix

Visualizing branch-cuts

Plot of Legendre function Q3

2z over a complex plane shows singularities, and

branch-cuts:

• ,

,

• LegendreQ3, 2, z 825 z215 z4

1z2

15 z 1 z2 ArcTanhz

16 Talk - Part I.nb

Q3

2z is continuous at branch line z 1 from below and at z 1 from above:

TableForm TableLimitLegendreQ3, 2, z, z p, Direction dir FullSimplify, dir, I, I, 1, p, 2, 2, TableHeadings "from above", "from below", "at the point", z 2, z 2 z 2 z 2 from above

148 3

45 Π 45 Log3

148 3

45 Π 45 Log3 from below

148 3

45 Π 45 Log3

148 3

45 Π 45 Log3 at the point

148 3

45 Π 45 Log3

148 3

45 Π 45 Log3

Talk - Part I.nb 17

Special functions in Mathematica

Numeric evaluation: multiple algorithms

Use of asymptotic expansions

Direct summation of the defining series vs. asymptotic expansion: BesselJ3.3, 10. ^ 3 0.00691901 WithΝ 3.3, x 1000., ListLogPlotTable 1 ^ k x 2Ν2 k k GammaΝ 1 k , k, 0, 1500

200 400 600 800 1000 1200 1400 1039 1056 10151 10246 10341

WithΝ 3.3, x 1000., MaxAbsTable 1 ^ k x 2Ν2 k k GammaΝ 1 k , k, 0, 1500 Log10 430.795020822201 WithΝ 33 10, x 1000``450, Sum 1 ^ k x 2Ν2 k k GammaΝ 1 k , k, 0, 1500 0.00691901025297857

18 Talk - Part I.nb

Accuracy 448.402

Divide-and-conquer techniques

Index doubling: AbsoluteTiming Fibonacci10 ^ 7 IntegerLength 0.1092002, 2 089 877 ChebyshevT10 ^ 7, 2 IntegerLength AbsoluteTiming 0.2028004, 5 719 476 v1 ChebyshevT10 ^ 4, 2; AbsoluteTiming 0., Null poly ChebyshevT10 ^ 4, x; AbsoluteTiming 0.7800014, Null v1 poly . x 2 True Use of divide-and-conquer summation at rationals (see [Haible, Papanikolaou], also [Karatsuba]) HermiteH10 ^ 5, 2 IntegerLength AbsoluteTiming 0.3276006, 243 338 LegendreP10 ^ 5, 3 IntegerLength AbsoluteTiming 0.4836008, 76 553

Truncated series for exponential: expn1

k0 n nk nk 1 k

ra HypergeometricPFQ10 000, 10 000, 1; AbsoluteTiming 0.0312001, Null NExp1, 35 700 ra 3.51338283258400315597019444219345045 1035 664 1.0 k . k 10 001 3.513031564557862 1035 664

Talk - Part I.nb 19

Numeric vs. symbolic

Numeric evaluation is used: HypergeometricPFQ1, 1, 1, 3 2, 2, 1`50 2.467401100272339654708622749969037783828424851810 Symbolic evaluation is used first, then use N: HypergeometricPFQ1, 1, 1, 3 2, 2, 1 Π2 4 N, 50 2.4674011002723396547086227499690377838284248518102

20 Talk - Part I.nb

Special functions in Mathematica

differentiation (see D, Derivative)

Expression for the derivative may be not unique, nicer one is chosen as a default: DBesselJΝ, x, x 1 2 BesselJ1 Ν, x BesselJ1 Ν, x Ν BesselJΝ, x x BesselJ1 Ν, x BesselJ1 Ν, x Ν BesselJΝ, x x FullSimplify True Special values of derivative: DBesselJΝ, x, Ν . Ν 0 BesselJ1,00, x

Comment

Alternative input:

Derivative1, 0BesselJ0, x BesselJ1,00, x FunctionExpand

Talk - Part I.nb 21

Special functions in Mathematica

series expansion (see Series also Limit, Residue)

Expansion at generic points

Genererating function of Catalan numbers: Series 2 1 1 4 x , x, 0, 12 1 x 2 x2 5 x3 14 x4 42 x5 132 x6 429 x7 1430 x8 4862 x9 16 796 x10 58 786 x11 208 012 x12 Ox13 Expansion of Aix at a generic point: SeriesAiryAix, x, a, 4 AiryAia AiryAiPrimea x a 1 2 a AiryAia x a2 1 6 AiryAia a AiryAiPrimea x a3 1 24 a2 AiryAia 2 AiryAiPrimea x a4 Ox a5

Expansion at poles

Explicit point: SeriesCotΠ x, x, 1, 4 1 Π x 1 1 3 Π x 1 1 45 Π3 x 13 Ox 15 Expansion of Euler -function at symbolic negative integer point: SeriesGammas, s, k, 1, Assumptions k Integers && k 0 1k k s k 1k PolyGamma0, 1 k k

• 1

6 k 1k Π2 3 PolyGamma0, 1 k2 3 PolyGamma1, 1 k s k Os k2

22 Talk - Part I.nb

Expansion at regular singular points

Wronskian of independent solutions of Bessel’s ODE: SeriesDet BesselJΝ, x BesselYΝ, x DBesselJΝ, x, x DBesselYΝ, x, x , x, 0, 2 x2 Ν 222 Ν CosΠ 1 Ν Gamma1 Ν x Π Gamma1 Ν Ox3 x2 Ν 222 Ν CosΠ Ν GammaΝ x Π Gamma2 Ν Ox3

• 2

Π x 1 4 Π 1 Ν Ν 1 4 Π Ν 1 Ν Gamma1 Ν 4 Π Gamma1 Ν GammaΝ 4 Π Gamma2 Ν x Ox2

• x2 Ν

22 Ν CosΠ Ν GammaΝ Π GammaΝ x 222 Ν 1 2 Ν CosΠ Ν GammaΝ x Π Ν 1 Ν GammaΝ Ox2

• x2 Ν 22 Ν CosΠ 1 Ν Gamma1 Ν Π Gamma1 Ν x

222 Ν 1 2 Ν CosΠ 1 Ν Gamma1 Ν x Π Ν 1 Ν Gamma1 Ν Ox2 FullSimplify x2 Ν Ox2 2 Π x Ox2

Compare:

WronskianBesselJΝ, x, BesselYΝ, x, x 2 Π x

Talk - Part I.nb 23

Expansion at branch point

Point z 1 is a regular singular point and a branch-point of Gauss’s hypergeo- metric function: SeriesHypergeometric2F1a, b, a b 1, z, z, 1, 1 Floor Arg1 z 2 Π

• 2 Π Gamma1 a b z 1

Gammaa Gammab Oz 12

• Gamma1 a b

Gamma1 a Gamma1 b 1 Gammaa Gammab Gamma1 a b 1 2 EulerGamma Π Log1 z PolyGamma0, 1 a PolyGamma0, 1 b z 1 Oz 12 Branch-cut corrections are artificial, but needed: SeriesHypergeometric2F1a, b, a b 1, 1 w, w, 0, 1 Gamma1 a b Gamma1 a Gamma1 b 1 Gammaa Gammab Gamma1 a b 1 2 EulerGamma Logw PolyGamma0, 1 a PolyGamma0, 1 b w Ow2

24 Talk - Part I.nb

Expansion at irregular singular points

SeriesBesselIΝ, x, x, Infinity, 3 Expand x

1 x

2 Π

• 1 4 Ν2 1

x 32

8 2 Π

• 9 40 Ν2 16 Ν4 1

x 52

128 2 Π O 1 x

72

• x

Π Ν

1 x

2 Π

• Π Ν 1 4 Ν2 1

x 32

8 2 Π

• Π Ν 9 40 Ν2 16 Ν4 1

x

52

128 2 Π O 1 x

72

Series Gammak Α Gammak Β , k, , 1 1 k

ΑΒ

1 Α Α2 Β Β2 2 k O 1 k

2

For further reference: Stirling’s formula

Generalized Stirling’s formula:

logx a x a 1 2 logx x 1 2 log2 Π

m1 p

1m1 m m 1 Bm1a xm xp

Withp 10, SeriesLogGammax a x a 1 2 Logx x 1 2 Log2 Π

• m1

p

1m1 m m 1 BernoulliBm 1, a 1 xm , x, , p Simplify True

Talk - Part I.nb 25

Special functions in Mathematica

functional expansion: to simpler, less general (see FunctionExpand)

Expand Meijer’s G-function MeijerG0, 0, , 1, 1, , x MeijerG0, 0, , 1, 1, , x FunctionExpand x 2 2 x 1 x Logx 1 x3 Expansion hypergeometric function in terms on simpler special functions: HypergeometricPFQ3 2, 2, 3, x HypergeometricPFQ 3 2 , 2, 3, x FunctionExpand 4 BesselI1, x

2

x Expansion of parabolic cylinder function in terms of elementary functions: ParabolicCylinderD3, x ParabolicCylinderD3, x FunctionExpand x2

4 6

2 x 2 2 x3 2 2 Special values: Assumingn Integers && n 1, FunctionExpandZeta2 n Refine 1n 212 n Π2 n BernoulliB2 n Gamma1 2 n

26 Talk - Part I.nb

Use FunctionExpand to expand holonomic DifferenceRoot/DifferentialRoot in special function: FunctionExpandDifferentialRootFunctiony, x, y ''x x yx, y0 c1, y '0 c2z c2 AiryAiz 323 c2 AiryBiz 3 316 Gamma 2 3 1 2 323 c1 AiryAiz Gamma 2 3 1 2 316 c1 AiryBiz Gamma 2 3

Talk - Part I.nb 27

Special functions in Mathematica

simplification (see FullSimplify)

Simplification blindly applies applicable transformations to subexpressions, keeping only those that reduce complexity. Identities with related functions

SinΠ Ν BesselYΝ, x BesselJΝ, x FullSimplify BesselJΝ, x CosΠ Ν ExpIntegralEiLogx FullSimplify LogIntegralx Gammaa Gamma1 b Hypergeometric1F1a, b, z Gamma1 a b z1b Gammab 1 Hypergeometric1F11 a b, 2 b, z FullSimplify Gammaa HypergeometricUa, b, z

Recurrence equations

Differentiation rule for the fraction of Bessel functions: D BesselJΝ 1, x BesselJΝ, x , x BesselJ1 Ν, x BesselJ1 Ν, x BesselJ1 Ν, x 2 BesselJΝ, x2 BesselJΝ, x BesselJ2 Ν, x 2 BesselJΝ, x Simplify`SimplifyRecurrence 1 1 2 Ν BesselJ1 Ν, x x BesselJΝ, x BesselJ1 Ν, x2 BesselJΝ, x2

Thus we established:

28 Talk - Part I.nb

g 'x 1 2 Ν 1 x gx gx2 . g Functionx, BesselJΝ 1, x BesselJΝ, x FullSimplify True Verify solutions of differential equations: x 1 x y ''x Γ Α Β 1 x y 'x Α Β yx . y Functionx, c1 Hypergeometric2F1Α, Β, Γ, x x1Γ c2 Hypergeometric2F11 Α Γ, 1 Β Γ, 2 Γ, x FullSimplify

Word of caution

FullSimplify may reject transformations that increase complexity of intermedi- ate expressions: PolyLog2, z PolyLog2, 1 z FullSimplify, 0 z 1 & PolyLog2, 1 z PolyLog2, z

Applying FunctionExpand forces those transformations to be applied:

PolyLog2, z PolyLog2, 1 z FunctionExpand FullSimplify, 0 z 1 &

Talk - Part I.nb 29

Special functions in Mathematica

support in solvers (Integrate, Sum, DSolve, etc.)

Used as output

Integrals: IntegrateExpz t t12 1 t12, t, 0, , Assumptions z 0 1 2 Π HypergeometricU 3 2 , 2, z FunctionExpand 1 2 z2 BesselK0, z 2 1 2 z2 BesselK1, z 2 Integrateta1 1 tca1 1 t zb, t, 0, 1, Assumptions z 1 && c a 0 Gammaa Gammaa c Hypergeometric2F1Regularizeda, b, c, z FunctionExpand Gammaa Gammaa c Hypergeometric2F1a, b, c, z Gammac Sums: SumBinomialn, k2 zk, k, 0, n 1 zn LegendrePn, 1 z 1 z SumBinomialn, k4 zk, k, 0, n HypergeometricPFQn, n, n, n, 1, 1, 1, z

30 Talk - Part I.nb

Differential equations: DSolvey ''x 3 x y 'x 4 x ^ 2 yx 0, y, x y Functionx, 2 x2 C1 HermiteH 4 5 , 5 2 x 2 x2 C2 Hypergeometric1F1 2 5 , 1 2 , 5 x2 2

• Recurrence equations:

RSolve2 n 1 fn n 1 fn 1 f1 1, f, n f Functionn, 3 22n Pochhammer1, 1 n Pochhammer 1

2 , 1 n

• Used as input

In integrals: IntegrateArcTanx BesselK0, z x, x, 0, , Assumptions z 0 1 8 MeijerG0, 1 2 , 1 2 , 0, 0, 0, , z2 4 N . z 2 0.211963 NIntegrateArcTanx BesselK0, 2 x, x, 0, 0.211963 Sums: SumCatalanNumberk CatalanNumbern 1 k, k, 0, n 1 CatalanNumbern Solving equations: ReduceBesselJ0, x BesselJ2, x && 0 x 2, x

Reduce::nint : Warning: Reduce used numeric integration to show that the solution set found is complete.

x RootBesselJ0, 1 BesselJ2, 1 &, 1.84118378134065930264362951364

Talk - Part I.nb 31

PlotBesselJ0, x BesselJ2, x, x, 0, 2

0.5 1.0 1.5 2.0 0.2 0.4 0.6 0.8 1.0

Solving differential equations: DSolvey ''x yx SinIntegralx, y, x y Functionx, C1 Cosx C2 Sinx 1 2 CosIntegral2 x Sinx Logx Sinx 2 Cosx2 SinIntegralx 2 Sinx2 SinIntegralx Cosx SinIntegral2 x

Neat example

Holonomic representation of a hypergeometric function: DifferenceRootReduce Hypergeometric2F1n 1 2, n 1 2, 2 n 1, z, n DifferenceRootFunctiony . ., n . ., 64 1 n . . 2 n . . y . .n . . 32 1 n . . 2 n . . 2 z y . .1 n . . 3 2 n . .2 z2 y . .2 n . . 0, y . .0 2 EllipticKz Π , y . .1 Hypergeometric2F1 3 2 , 3 2 , 3, zn

Notice that for z 2 the middle term of the recurrence relation vanishes. Use FunctionExpand to try to solve the obtained recurrence equation:

32 Talk - Part I.nb

FunctionExpand . z 2

• n Π

2

EllipticK2 Gamma 3 4

2

Gamma1 n Π Gamma 3 4 n 2

2

• 1n

n Π 2

EllipticK2 Gamma 3 4

2

Gamma1 n Π Gamma 3 4 n 2

2

• 4

n Π 2

EllipticE2 Gamma 5 4

2

Gamma1 n Π Gamma 3 4 n 2

2

• 4 1n

n Π 2

EllipticE2 Gamma 5 4

2

Gamma1 n Π Gamma 3 4 n 2

2

Hypergeometric2F1n 1 2, n 1 2, 2 n 1, 2 FullSimplify, n Integers && n 0 Hypergeometric2F1 1 2 n, 1 2 n, 1 2 n, 2 2 1n n Cos n Π 2 EllipticK2 Gamma 3 4

2

4 EllipticE 2 Gamma 5 4

2

Sin n Π 2

• Π Gamma 3

4 n 2

2

Verify:

Table, n, 0, 6 FullSimplify True, True, True, True, True, True, True

Talk - Part I.nb 33

Defining your own function

Define a function as symbolic expression in terms of known functions

Wave function of harmonic oscillator Ψn_, x_ : x2

2 HermiteHn, x

2n Π12 n

Comment

Advantage: solvers will see the input in terms of objects they know how to deal with.

IntegrateΨ4, x x ^ 2 Ψ4, x, x, , 9 2

Define atomic special function Let f x x cosxsinx

x2

.

Clearf; f0 : 0 fPi k_. ; IntegerQk : 1 ^ k Pi k

Evaluation for reals and complexes, using the form with no division by zero problem:

fx_Real : x 3 HypergeometricPFQ, 5 2, x2 4 fz_Complex : z 3 HypergeometricPFQ, 5 2, z2 4 f : Derivative1f : Functionz, 2 z Cosz 2 z2 Sinz z3

• Comment

Disadvantage of such an approach is that solvers see this function as a black box:

Integratefx, x fx x Integratex Cosx Sinx x ^ 2, x Sinx x

34 Talk - Part I.nb

Thank you

Initialization

Talk - Part I.nb 35