Special Functions in Mathematica

Oleksandr Pavlyk

Wolfram Research Inc.

Abstract: Special functions are often at heart of computations with Mathematica. We outline evolutionary improvements of legacy special functions, and present special functions new in upcoming Mathematica 6.0. Their usage will be exemplified through applications. New functions covered will include Kelvin functions, Whittaker functions, parabolic cylinder D, spherical Bessel functions, functions of spheroidal family, Bell numbers and Bell polynomials.

Initialization

SetDirectory[DirectoryName[ToFileName["FileName"/.NotebookInformation[EvaluationNotebook[]]]]];

<<"MeijerGFormDos.m"

<<"FractionalDerivative.m"

ClearAll[MyForm];

MyForm /: MakeBoxes[MyForm[t_Times], form_]:=Module[{t2,t3,t4,var},
t2 = (List @@ t);
t2 = {Times @@ Select[t2,NumericQ], Times @@ Complement[t2,Select[t2,NumericQ]]} ;
t3 =  (Unique[var, Temporary]& /@ t2);
t4 = Evaluate[(Thread[(ToString /@ t3)->(ToBoxes[#,form]& /@ t2)])];
ToBoxes[(Times @@ t3), form]/.t4
]

MyForm /: MakeBoxes[MyForm[t_], form_]:=MakeBoxes[t, form]

MeijerG; SetAttributes[MeijerGTools`MeijerGToSum,HoldAll];

MeijerGTools`MeijerGToSum[MeijerG[w___],k_]:=(MeijerG; System`MeijerGDump`MeijerGToSums[MeijerG[w],k])

VectorLength[x_List]:=Sqrt[Total[x^2]]

PlotSinc[Interval[{p1_,p2_}]]:=Module[{t},
t= N[Sinc[Interval[{p1,p2}]],30]/.Interval->Identity;
Labeled[Show[{
Plot[Evaluate@Join[{Sinc[x]},t],{x,-10,10},PlotStyle->{{Black},{Blue},{Purple}}],
Plot[Sinc[x],{x,p1,p2},PlotStyle->{Thickness[.005],Red}]
}],Grid[{{HoldForm[Sinc[Interval[{p1,p2}]]]->Sinc[Interval[{p1,p2}]]},{HoldForm[N@Sinc[Interval[{p1,p2}]]]->N@Sinc[Interval[{p1,p2}]]}}],Right]
]

Oleg Igorevich Marichev

Collaboration

• Improvements to Hypergeometric and MeijerG evaluations (logarithmic cases)

N[MeijerG[ {{},{}},{{1,1,1},{}},1/GoldenRatio],200]

MeijerGTools`MeijerGToSum[MeijerG[ {{},{}},{{1,1,1},{}}, z,1],k]

• Ongoing work in Hypergeometric expansion

• Continual discussions

Idea generation

• PowerExpand support of assumptions

{(Sqrt[x^2]),PowerExpand[ (Sqrt[x^2]),Assumptions->x<0], PowerExpand[(Sqrt[(x^2)])]}

%/.x->-1

• Series on branch-cuts

Series[ ArcCos[x],{x,2,4}]

Series[ ArcCos[x],{x,2,4},Assumptions->(x-2)>0]

Series[ ArcCos[x],{x,2,4},Assumptions->(0<Arg[x-2]<Pi/2)]

• Fractional differentiation:

ℱD[ 1, {x, p}]

Table[ ℱD[ 1, {x, p}],{p,-2,2,1/2}]

ℱD[ ArcSin[x],{x,p}]

Outline of the talk

Spheroidal functions

{SpheroidalPS[n,m,γ,x], SpheroidalQS[n,m,γ,x], SpheroidalS1[n,m,γ,x], SpheroidalS2[n,m,γ,x], SpheroidalEigenvalue[n,m,γ]}

Bessel family functions

{HankelH1[ν,x], HankelH2[ν,x], KelvinBer[ν,x], KelvinBei[ν,x], KelvinKer[ν,x], KelvinKei[ν,x], SphericalBesselJ[ν,x], SphericalBesselY[ν,x],SphericalHankelH1[ν,x],SphericalHankelH2[ν,x]}

Confluent hypergeometric family functions

{WhittakerM[k,m,x], WhittakerW[k,m,x], ParabolicCylinderD[ν,x]}

Number theory functions

{ BellB[n], BellB[n,x], BernoulliB[n] }

Sinc

FunctionExpand[Sinc[x],x!=0]

Wavelets:

Integrate[ ((Sin[Pi(x-y)])/((x-y)Pi))((Sin[Pi(z-y)])/((z-y)Pi)),{y,-∞,∞}, GenerateConditions->False]

Plot[Sinc[x],{x,-10 Pi,10 Pi},PlotRange->All]

{Sinc[Pi/7],FunctionExpand[Sinc[Pi/7]]}

Sinc and Interval

PlotSinc[Interval[{1,3}]]

PlotSinc[Interval[{-7,3}]]

Spheroidal special functions

These are solutions arising from separation of variables in Helmholtz differential equation in spheroidal coordinates

Spheroidal coordinates

Prolate

{x,y,z}->{ (Sqrt[(1-(η^2))((ζ^2)-1)])Cos[ϕ] ,(Sqrt[(1-(η^2))((ζ^2)-1)])Sin[ϕ], η ζ}

ParametricPlot[ {η ζ,(Sqrt[(1-η^2)(ζ^2-1)])}, {η,-1,1},{ζ,1,10},AxesLabel->{"z","r"}]

Oblate

{x,y,z}->{ (Sqrt[(1-(η^2))((ζ^2)+1)])Cos[ϕ] ,(Sqrt[(1-(η^2))((ζ^2)+1)])Sin[ϕ], η ζ}

ParametricPlot[ {η ζ,(Sqrt[(1-η^2)(ζ^2+1)])}, {η,-1,1},{ζ,0,10},AxesLabel->{"z","r"}]

Spheroidal special functions

Spheroidal differential equations for radial functions

(∂_ζ) ((1-(ζ^2))f'[ζ]) + ( λ + (γ^2)(1-(ζ^2))- ((m^2)/(1-(ζ^2))))f[ζ]==0

Solutions  f[ζ]==(Subscript[c, 1])SpheroidalS1[n,m,γ,x]+(Subscript[c, 2]) SpheroidalS2[n,m,γ,x]

Radial functions plots

Labeled[Plot[ Evaluate@Table[ Tooltip[SpheroidalS1[n,0,1,x],"n="<>ToString[n]],{n,0,4}],{x,-10,10},PlotRange->All],"Radial spheroidal functions "<>ToString[SpheroidalS1[n,0,1,x],StandardForm],Top]

Labeled[Plot[ Evaluate@Table[ Tooltip[SpheroidalS2[n,0,1,x],"n="<>ToString[n]],{n,0,4}],{x,1+10^-3,2}],"Radial spheroidal functions "<>ToString[SpheroidalS2[n,0,1,x],StandardForm],Top]

Spheroidal differential equations for angular functions

(∂_η) ((1-(η^2))f'[η]) + ( λ + (γ^2)(1-(η^2))- ((m^2)/(1-(η^2))))f[η]==0

Solutions  f[ζ]==(Subscript[c, 1])SpheroidalPS[n,m,γ,x]+(Subscript[c, 2]) SpheroidalQS[n,m,γ,x]

Angular functions plots

Labeled[ Plot[ Evaluate@Table[ Tooltip[SpheroidalPS[n,0,1,x],"n="<>ToString[n]],{n,0,6}],{x,-1,1}],"Angular spheroidal functions "<>ToString[SpheroidalPS[n,0,1,x],StandardForm],Top]

Labeled[ Plot[ Evaluate@Table[ Tooltip[SpheroidalQS[n,0,1,x],"n="<>ToString[n]],{n,0,3}],{x,-1+10^-3,1-10^-3}],"Angular spheroidal functions "<>ToString[SpheroidalQS[n,0,1,x],StandardForm],Top]

Spheroidal special functions

Relation to spherical functions (angular)

Series[ SpheroidalPS[n,m,γ,x],{γ,0,3}]

Series[ SpheroidalQS[n,m,γ,x],{γ,0,3}]

Relation to spherical functions (radial)

Normalization is fixed by asymptotics at γ z -> ∞:  SpheroidaS1[n,m,γ,z] -> (1/(γ z))Sin[γ z-n π/2], SpheroidaS2[n,m,γ,z] -> -(1/(γ z))Cos[γ z-n π/2]

1- Divide @@ With[{γ=1.`30,z=10^9,n=3,m=1},{SpheroidalS1[n,m,γ,z],(1/(γ z))Sin[γ z-(π/2)n]}]

1-Divide @@ With[{γ=1.`30,z=10^9,n=3,m=1},{SpheroidalS2[n,m,γ,z],-(1/(γ z))Cos[γ z-(π/2)n]}]

Special cases

for integer n spheroidal functions become elementary if γ=n π/2 and m=+1

Manipulate[StyleForm[MyForm@SpheroidalPS[n,1,((n π)/2),x],FontSize->18],{n,Range[10]},ControlType->Slider]

Spheroidal special functions

Prolate angular spheroidal functions SpheroidalPS[n,0,γ,x] are eigenfunctions of self-adjoint operator A:

A[SpheroidalPS][x]==(γ/π)(∫_(-1)^(-1))Sinc[γ (x-y)]SpheroidalPS[n,0,γ,y]dy == (Subscript[λ, n]) SpheroidalPS[n,0,γ,x]

Labeled[
Plot[
Block[{γ=1,n=0},γ/Pi NIntegrate[Sinc[γ (x-y)] SpheroidalPS[n,0,γ,y],{y,-1+(10^(-30)),1},WorkingPrecision->30,PrecisionGoal->15]/SpheroidalPS[n,0,γ,x]],{x,-1,1},WorkingPrecision->15,MaxRecursion->2,PlotPoints->4],StyleForm["A[SpheroidalPS]/SpheroidalPS for n=0, γ=1", "DisplayFormula"],{Left}]

Clear[λ];λ[n_Integer, γ_]:=(λ[n,γ]=Block[{x=1/5,r},r=γ/Pi NIntegrate[Sinc[γ (x-y)] SpheroidalPS[n,0,γ,y],{y,-1+(10^(-30)),1},WorkingPrecision->30,PrecisionGoal->20]/SpheroidalPS[n,0,γ,x];N[r,20]])

λ[#,1]& /@ Range[0,6]

ListPlot[{#,λ[0,#]}& /@ Range[1/10,1,1/10],PlotJoined->True,Mesh->All,AxesLabel->{γ,(Subscript[λ, 0])}]

Bessel family functions

Manipulate[Labeled[Plot[f[1/3,x],{x,(10^(-3)),10}],f[1/3,x]],{f,{SphericalBesselJ, SphericalBesselY, KelvinBer, KelvinBei, KelvinKer, KelvinKei}},ControlType->Checkbox]

Wire resistance and skin effect

Wire resistance increases with frequency due to skin effect. For circular cross-section wires resistance per unit length as a function of skin depth:

Res[q_]:=q/2 ( KelvinBer[0,q]KelvinBei'[q]-KelvinBei[0,q]KelvinBer'[q])/(KelvinBer'[q]^2+KelvinBei'[q]^2)

Plot[Res[q],{q,0.001,4},PlotRange->{.8,1.8}, ImageSize->300]

Series[Res[q],{q,0,20}]

Light diffraction

Light diffraction by an infinite conducting cylinder of unit radius:

ℰ[k,r,ϕ]==(e^(-i k r Cos[ϕ]))-((BesselJ[0,k R] HankelH2[0,k r])/(HankelH2[0,k R]))-2 (∑_(k=1)^∞)(((((-i))^m) BesselJ[m,k R] HankelH2[m,k r])/(HankelH2[m,k R]))

c[m_,kr_]:=(c[m,kr]=(-I)^m (BesselJ[m,kr]/HankelH2[m,kr]));
h[m_, kr_]:=(h[m, kr]=HankelH2[m,kr]);
ℰ[k_,r_,ϕ_,n_]:=Exp[-I k r Cos[ϕ]]+With[{R=1.},-c[0, k R]h[0, k r]-2 Sum[c[m, k R]h[m, k r] Cos[m ϕ],{m,1,n}]]

ListDensityPlot[Flatten[Outer[ {#1 Cos[#2], #1 Sin[#2], Abs[2 ℰ[2 Pi,#1,#2,12]]}&,Range[1,5,0.1],Range[0,2 Pi,0.06]],1],RegionFunction->(#1^2+#2^2>1&)]

Clear[c,h,ℰ];

New functions of confluent hypergeometric family

FunctionExpand[#[n,m,x]]& /@ {WhittakerM, WhittakerW}

Plot[{Tooltip[WhittakerM[2,1,x]],Tooltip@WhittakerW[2,1,x]},{x,1,5}]

Differential equation

Simplify[((∂_(x,x))#1+#1 (-(((m^2)-(1/4))/(x^2))+(k/x)-(1/4))&)[WhittakerW[k,m,x]]]

%/.(#->RandomReal[{-5,5},WorkingPrecision->100]& /@ {k,m,x})

FunctionExpand[
ParabolicCylinderD[n,x]==((1 )/(Sqrt[x]))(2^((n/2)+(1/4)))WhittakerW[(n/2)+(1/4),-(1/4),((x^2)/2)] 
]

Plot[Evaluate@(Tooltip/@ ParabolicCylinderD[Range[3],x]),{x,-5,5},PlotRange->All]

FunctionExpand[ParabolicCylinderD[-1/2,x]]

Manipulate[Tooltip[Simplify@FunctionExpand[ParabolicCylinderD[n,x]],n],{n,Range[5]},ControlType->Slider]

Green functions of Coulomb potential

Whittaker functions appear as Green functions in some non-relativistic QM problems in coordinate representation. In Coulomb problem

CoulombGreenFunction[r1_List,r2_List,ℰ_]:=Block[{ν,a1,a2,δ,𝓋ℓ=VectorLength},ν=1/(Sqrt[ℰ]);δ=𝓋ℓ[r2-r1];{a1,a2}= 𝓋ℓ[r2]+𝓋ℓ[r1]+{1,-1}δ;
-Gamma[1-I ν](1/(4 π δ))Det[(({{WhittakerW[I ν,(1/2),-((I a2)/ν)], WhittakerM[I ν,(1/2),-((I a1)/ν)]}, {Derivative[0,0,1][WhittakerW][I ν,(1/2),-((I a2)/ν)], Derivative[0,0,1][WhittakerM][I ν,(1/2),-((I a1)/ν)]}}))]
]

Verify numerically that it solve Schrödinger equation with Coulomb potential:

ReplaceAll[
((∂_({x1,2}))#1+(∂_({x2,2}))#1+(∂_({x3,2}))#1&)[CoulombGreenFunction[{x1,x2,x3},{y1,y2,y3},ℰ]]+((2/(VectorLength[{x1,x2,x3}]))+ℰ)CoulombGreenFunction[{x1,x2,x3},{y1,y2,y3},ℰ],
Append[(#->RandomReal[{1,20},WorkingPrecision->40]& /@ {x1,x2,x3,y1,y2,y3}),ℰ->-3/7]]

Poles correspond to bound state energies

Plot[Evaluate@CoulombGreenFunction[{1,0,0},{0,0,1/2},ℰ],{ℰ,-3/2,-1/100},WorkingPrecision->30]

Number theory functions

The number of ways a set of n elements can be partitioned into nonempty subsets is called a Bell number

Prime BellB numbers up to 60:

Position[BellB[Range[60]],x_/;PrimeQ[x]]

Corresponding BellB numbers

BellB @@@ %

Next prime Bell number (Subscript[B, 2841]) was proved  prime by by I. Larrosa Canestro in 2004 after 17 months of computation using the elliptic curve primality proving program PRIMO. (fact taken from MathWorld)

Timing[p = BellB[2841]; IntegerDigits[p]//Length]

p

Number theory functions

BernoulliB got improved for large numbers using algorithm of Bernd C. Kellner

(p = BernoulliB[5000];)//Timing

Length[IntegerDigits[Numerator[p]]]

7 24 3

BernoulliB[5 10^6] took about 500 GHz-hours of computation:

bn =BernoulliB[500000]; //AbsoluteTiming

{5191.0761519`11.166802393422284 Second,Null}

bn2 = BernoulliB[2000000]; //AbsoluteTiming

{118613.558123`12.525679327393227 Second,Null}

bn3=BernoulliB[5000000]; //AbsoluteTiming

{710866.8270827`13.303333241651519 Second,Null}

IntegerDigits[Numerator[bn3]]//Length

27332507

Row[{Button["Set SlideShow",SetOptions[EvaluationNotebook[],ScreenStyleEnvironment->"SlideShow"]],Button["Set Working",SetOptions[EvaluationNotebook[],ScreenStyleEnvironment->"Working"]]}]