Chapter 6
Special Functions and Orthogonal Polynomials

6.1 Hermite Polynomials

The Hermite polynomials can be obtained in a number of ways within Mathematica, e.g. by calling the in-built function HermiteH.  We can also solve the recursion relations (6.14) to generate the Hermite polynomials directly by using the ability of Mathematica to recognize and solve such equations.  We first construct (6.14), retaining both the lambda-dependence and the n-dependence in the expansion coefficients:

a[k_,n_]:=-2((k-n+2)/(n(n-1)))a[k,n-2]

The initial values for these equations are determined by  stipulating that a0 must 
vanish if k is an odd integer, and that a1 must vanish if k is an even integer. Furthermore, to obtain the standard form of the H_n in (6.22), whereby the coefficient of x^n is set equal to 2^n, the coefficients a0 and a1 must be assigned the values that are calculated in Problem 5:

a[k_,0]:=If[EvenQ[k],(-1)^(k/2) k!/(k/2)!,0]
a[k_,1]:=If[EvenQ[k],0,2(-1)^((k-1)/2) k!/((k-1)/2)!]

The Hermite polynomials are now defined by using the series solutions in (6.22).  Both series can be combined into a single expression by using conditional statements to distinguish between even-order and odd-order polynomials:

One[k_] :=If[EvenQ[k],0,1]
kMax[k_]:=If[EvenQ[k],k/2,(k-1)/2]
H[k_,x_]:=Sum[a[k,2j+One[k]]x^(2j+One[k]),{j,0,kMax[k]}]

Expressions for the Hermite polynomial of order k can now be generated by calling the function H[k,x]:

Do[Print[H[k,x]],{k,0,10}]

Figure 6.1

The graphs of the Weber-Hermite functions shown in Figure 6.1 are generated by first constructing the (normalized) functions using (6.24) and the normalization derived in Problem 39:

Psi[n_,x_]:=(Sqrt[(1/(Sqrt[Pi]2^n n!))]*
            Exp[-(1/2)x^2]HermiteH[n,x])

The individual panels of the figure are then obtained with Plot as follows:

Do[Show[{Plot[Psi[k,x],{x,-6,6},PlotRange->{-0.8,0.9},
     Ticks->{Automatic,Range[-0.8,0.8,0.4]},
     AxesLabel->{x,psi},DisplayFunction->Identity],
     Graphics[{Text["n=",{5,0.6}],Text[k,{5.5,0.6}]}]},
     DisplayFunction->$DisplayFunction],
     {k,0,9}]

6.2 Legendre Polynomials

Figure 6.2

The spherical polar plots of the first few Legendre polynomials in this figure are generated with ParametricPlot3D.  After loading this package, the coordinate axes are first defined in the graphics object CoordinatesAxes:

<<Graphics`ParametricPlot3D`

CoordinateAxes[t_]:=
     ParametricPlot3D[{{u,0,0},{0,u,0},{0,0,u}},
     {u,-t,t},DisplayFunction->Identity]

The spherical polar plot of the Legendre polynomial L[l,x] is then also defined as a graphics object with the aid of the in-built function LegendreP,

y[l_]:=LegendreP[l,Cos[theta]];
LegendrePolynomial[l_]:=
     ParametricPlot3D[{Abs[y[l]]Sin[theta]Cos[phi],
     Abs[y[l]]Sin[theta]Sin[phi],Abs[y[l]]Cos[theta]},
     {theta,0,Pi,Pi/30},{phi,0,2Pi,Pi/15},
     DisplayFunction->Identity]

and the coordinate axes and the Legendre polynomial are combined with Show:

LegendrePlot[l_,t_]:=Show[CoordinateAxes[t],
     LegendrePolynomial[l],Boxed->False,Axes->False,
     AmbientLight->GrayLevel[1],LightSources->{},
     PlotRange->All,DisplayFunction->$DisplayFunction]

Each panel of the figure can now be obtained by using LegendrePlot with the appropriate arguments:

LegendrePlot[0,1.25]

LegendrePlot[1,1.25]

LegendrePlot[2,1.25]

LegendrePlot[3,1.25]

6.3 Legendre Functions and Spherical Harmonics

Figure 6.3

The panels in this figure are obtained by following similar steps to those used to generate Figure 6.2.  The required packages are first loaded

<<Graphics`ParametricPlot3D`

and the coordinate axes are defined as before:

CoordinateAxes[t_]:=
     ParametricPlot3D[{{u,0,0},{0,u,0},{0,0,u}},
     {u,-t,t},DisplayFunction->Identity]

The remaining steps are now carried out as before:

SurfaceHarmonic[y_]:=
     ParametricPlot3D[{Abs[y]Sin[theta]Cos[phi],
     Abs[y]Sin[theta]Sin[phi],Abs[y]Cos[theta]},
     {theta,0,Pi,Pi/30},{phi,0,2Pi,Pi/15},
     DisplayFunction->Identity]

YlmPlot[y_,t_]:=Show[CoordinateAxes[t],
     SurfaceHarmonic[y],Boxed->False,Axes->False,
     AmbientLight->GrayLevel[1],LightSources->{},
     PlotRange->All,DisplayFunction->$DisplayFunction]

YlmPlot[Sin[theta]Cos[phi],1.5]

YlmPlot[Sin[theta]Sin[phi],1.5]

YlmPlot[(Sin[theta])^2 Cos[2 phi],1.5]

YlmPlot[(Sin[theta])^2 Sin[2 phi],1.5]

YlmPlot[Sin[theta]Cos[theta]Cos[phi],0.75]

YlmPlot[Sin[theta]Cos[theta]Sin[phi],0.75]

The alert reader will have noticed that the diagrams produced here differ by a rotation of Pi/2 in the x-y plane from those in the book.  This is simply a matter of relabelling the x and y axes (since absolute values are being plotted).

6.4 Bessel Functions

Although the Bessel functions of the first and second kinds are in-built functions in Mathematica (BesselJ and BesselY, respectively), we can use the series representations (6.86) and (6.98) to investigate their convergence.  Before doing so, we first examine the convergence of the limit in (6.100), which is Euler's constant.  We first construct the partial sum in (6.100) and then take the limit k->Infinity:

S[k_]:=Sum[1/n,{n,1,k}]

Do[Print[N[S[k]-Log[k],10]],{k,1,2000,200}]

The convergence is seen to be very slow, as even after 2000 terms, the sequence 
has converged only to the first three decimal places.  This result is to be compared
with the `exact' value obtained from EulerGamma to ten decimal places:

N[EulerGamma,10]

This slow convergence of our construction makes it much more convenient to 
use the in-built function PolyGamma to calculate the logarithmic derivative of the
Gamma function instead of performing explicitly the operations in (6.99).

To examine the convergence of the series representation for Y_n(x), we first construct the partial sums of (6.98):

Y[n_,x_,k_]:=(2/Pi)BesselJ[n,x]Log[x/2]-
     (1/Pi)Sum[((n-j-1)!/j!)(x/2)^(2j-n),{j,0,n-1}]-
     (1/Pi)Sum[((-1)^j/(j!(n+j)!))(PolyGamma[j+1]+
          PolyGamma[j+n+1])(x/2)^(2j+n),{j,0,k}]

This series is now evaluated as a function of the cutoff p in the summation of the last term on the right-hand side of (6.98)  for values of the argument of x=1 and x=5 for Y_2 and Y_8:

Do[Print[N[Y[2,1,i],10]," ",N[Y[8,1,i],10]],{i,1,5}]

Do[Print[N[Y[2,5,i],10]," ",N[Y[8,5,i],10]],{i,1,5}]

The convergence for Y_2(1) and Y_8(1) is seen to be very rapid.  This is because the first two terms on the right-hand side of (6.98) dominate, with the contribution from the series providing a relatively small correction.  However,with increasing n, the convergence is faster for larger arguments, as can be seen by comparing the sequence of partial sums Y_2(5) and Y_8(5).

6.5 Generating Functions

Special Functions from Generating Functions

Special functions and orthogonal polynomials can be obtained easily from their 
generating function in Mathematica.  For the Legendre polynomials, we first 
construct the right-hand side of (6.115):

G[x,t]=1/Sqrt[1+t^2-2 t x];

The Legendre polynomial of n-th order is obtained by differentiating this function n times with respect to t, dividing by n! and then setting t=0:

P[n_,x_]:=Together[D[G[x,t],{t,n}]/n!/.t->0];

The Legendre polynomials can now be generated directly from P[n,x].  For
example, to obtain P_5(x), we simply enter

P[5,x]

Recursion Relations

Mathematica can be used to verify the recursion relations explicitly.  Equation 
(6.119) is first entered by writing the k-th recursion relation in terms of the in-built function LegendreP as

Recursion[k_]:=Simplify[k LegendreP[k-1,x]+
     (k+1)LegendreP[k+1,x]-(2k+1) x LegendreP[k,x]==0]

We can now test this for any value of k:

Recursion[10]

The recurrence relation (6.119) can also be used to generate the hierarchy of
Legendre polynomials by specifying only the `initial conditions' P_0(x)=1 and P_1(x)=x.

P[0,x_]:=1

P[1,x]:=x

P[n_,x_]:=Together[((2n-1)x P[n-1,x]-(n-1)P[n-2,x])/n]

Any of the Legendre polynomials can now be generated directly from the 
recursion relations.  For example,

P[5,x]