Solutions to Differential Equations: Special Functions

Visualizing a few special functions

DSolve[
      x^2 y''[x]  +  
          x y'[x]  +  (x^2  -  ν^2) y[x]  == 0,  
      y[x], x]

MyPlotStyle[HowMany_Integer]  :=  
      Table[{Hue[0.66*i/HowMany], Thickness[0.01]}, {i, 0, 
          HowMany}]

Plot[{BesselJ[0, x],  BesselY[0, x]}, {x, 0, 20},  
      PlotStyle -> MyPlotStyle[2]]

Plot[{BesselJ[5, x],  BesselY[5, x]}, {x, 0, 20},  
      PlotStyle -> MyPlotStyle[2]]

Plot[{BesselJ[1/2, x],  BesselY[1/2, x]}, {x, 0, 20},  
      PlotStyle -> MyPlotStyle[2]]

DSolve[(1 - x^2 ) y''[x]  -  
          2  x y'[x]  +  (n (n + 1)) y[x]  == 0,  
      y[x], x]

DSolve[(1 - x^2 ) y''[x]  -  
          2  x y'[x]  +  (n (n + 1)  -  m^(2/(1 - x^2))) y[
              x]  == 0,  y[x], x]

Plot[{LegendreP[0, x],  LegendreQ[0, x]}, {x, -1, 1},  
      PlotStyle -> MyPlotStyle[2]]

Plot[{LegendreP[1, x],  LegendreQ[1, x]}, {x, -1, 1},  
      PlotStyle -> MyPlotStyle[2]]

Plot[{LegendreP[2, x],  LegendreQ[2, x]}, {x, -1, 1},  
      PlotStyle -> MyPlotStyle[2]]

Plot[Evaluate[Table[LegendreP[i, x], {i, 0, 10}]], {x, -1, 1},  
      PlotStyle -> MyPlotStyle[11]]

Plot[Evaluate[Table[LegendreQ[i, x], {i, 0, 10}]], {x, -1, 1},  
      PlotStyle -> MyPlotStyle[11]]

Using Special Functions to visualize the the eigensolutions for an electron bound to a fixed proton--The Hydrogen Atom

DSolve[
      x y''[x]  +  (q   +  1  -  x) y'[x]  +  
          p  y[x]  == 0,  y[x], x]

Plot[LaguerreL[4, 1, x], {x, -5, 15}]

Eigenfunctions for Hydrogen:

Preliminary Definitions:
n will be the first quantum number n =1,2,3
L will be the second quantum number L=0,1,2,... (n-1)
m will be the third quantum number m = -L,-L+1,... 0,1,2, +L
These rules come from restrictions on the special functions that make up the eigenfucntions for the H-atom

Z = 1;  a0  =  1;
κ[
        n_]  :=  (( Z)/(a0  n))
ρ[n_,  r_]  :=  2 κ[n] r

    
A[n_ ,  
        L_]  :=  (Sqrt[Factorial[n - L  -  1]/(2 n (((Factorial[n  +  L]))^3))])
R[n_ ,  L_,  
        r_]  :=  (((2 κ[n]))^((3/(2 )))) A[n, 
          L] ρ[n, r] Exp[((-ρ[n, r])/2)] LaguerreL[n - L - 1, 
          2 L  +  1, ρ[n, r]]

The eigenfunctions for the Hydrogen Atom can now be written down.  The eigenfunctions will relate to the above definitions as well as the Spherical Harmonics are related to Legendre functions and its derivatives.

HAtom[n_ ,  L_ ,  m_ , {r_ ,  θ_,  ϕ_}]  :=  
      R[n, L, r] SphericalHarmonicY[L, m, θ, ϕ]

<< Calculus`VectorAnalysis`

CoordinatesFromCartesian[{x, y, z}, Spherical]

HAtom[2, 1, 1, 
      CoordinatesFromCartesian[{1, 2, 1}, Spherical]]

N[HAtom[2, 1, 1, 
        CoordinatesFromCartesian[{1, 2, 1}, Spherical]]]

<< Graphics`ContourPlot3D`

Abs[
        HAtom[2, 1, 1, 
          CoordinatesFromCartesian[{0, 1,  .5}, Spherical]]]^2

Plot a surface of constant probability density. Only an octant can be plotted, the numerical methods do not behave well on the coordinate "zero value" planes.

p1 = 
      ContourPlot3D[
        Abs[HAtom[2, 1, 1, 
                CoordinatesFromCartesian[{x, y, z}, Spherical]]]^2 - 
          0.00005082450653152705, {x, 0.01, 5}, {y, 0.01, 5}, {z, 0.01, 
          5}]

p1 = 
      ContourPlot3D[
        Abs[HAtom[2, 1, 1, 
                CoordinatesFromCartesian[{x, y, z}, Spherical]]]^2 - 
          0.00005082450653152705, {x, 0.01, 5}, {y, 0.01, 5}, {z, 0.01, 5}, 
        MaxRecursion -> 2]

p2 = 
      ContourPlot3D[
        Abs[HAtom[2, 1, 1, 
                CoordinatesFromCartesian[{x, y, z}, Spherical]]]^2 - 
          0.00005082450653152705, {x, 0.0001, 5}, {y, 0.0001, 
          5}, {z, -0.0001, -5}, MaxRecursion -> 2]

However, we can always put the octants back together like a puzzle

Show[p1, p2]

Write a function to place all the octants on one graph

HAtomShow[n_ ,  L_,  m_,  probdensity_]  :=  
      Module[{nwu, nwd, swu, swd, neu, ned, seu, 
          sed}, 
nwu  =  
          ContourPlot3D[
            Abs[HAtom[n, L, m, 
                    CoordinatesFromCartesian[{x, y, z}, Spherical]]]^2 - 
              probdensity, {x, -0.0001, -5}, {y, 0.0001, 5}, {z, 
              0.0001, 5}, MaxRecursion -> 2, 
            DisplayFunction -> Identity]; 
nwd =  
          ContourPlot3D[
            Abs[HAtom[n, L, m, 
                    CoordinatesFromCartesian[{x, y, z}, Spherical]]]^2 - 
              probdensity, {x, -0.0001, -5}, {y, 0.0001, 
              5}, {z, -0.0001, -5}, MaxRecursion -> 2, 
            DisplayFunction -> Identity]; 
swu  =  
          ContourPlot3D[
            Abs[HAtom[n, L, m, 
                    CoordinatesFromCartesian[{x, y, z}, Spherical]]]^2 - 
              probdensity, {x, -0.0001, -5}, {y, -0.0001, -5}, {z, 0.0001, 5}, MaxRecursion -> 2, 
            DisplayFunction -> Identity]; 
swd =  
          ContourPlot3D[
            Abs[HAtom[n, L, m, 
                    CoordinatesFromCartesian[{x, y, z}, Spherical]]]^2 - 
              probdensity, {x, -0.0001, -5}, {y, -0.0001, -5}, {z, -0.0001, -5}, MaxRecursion -> 2, 
            DisplayFunction -> Identity]; 
neu  =  
          ContourPlot3D[
            Abs[HAtom[n, L, m, 
                    CoordinatesFromCartesian[{x, y, z}, Spherical]]]^2 - 
              probdensity, {x, 0.0001, 5}, {y, 0.0001, 5}, {z, 0.0001, 5}, 
            MaxRecursion -> 2, 
            DisplayFunction -> Identity]; 
ned =  
          ContourPlot3D[
            Abs[HAtom[n, L, m, 
                    CoordinatesFromCartesian[{x, y, z}, Spherical]]]^2 - 
              probdensity, {x, 0.0001, 5}, {y, 0.0001, 
              5}, {z, -0.0001, -5}, MaxRecursion -> 2, 
            DisplayFunction -> Identity]; 
seu  =  
          ContourPlot3D[
            Abs[HAtom[n, L, m, 
                    CoordinatesFromCartesian[{x, y, z}, Spherical]]]^2 - 
              probdensity, {x, 0.0001, 5}, {y, -0.0001, -5}, {z, 
              0.0001, 5}, MaxRecursion -> 2, 
            DisplayFunction -> Identity]; 
sed =  
          ContourPlot3D[
            Abs[HAtom[n, L, m, 
                    CoordinatesFromCartesian[{x, y, z}, Spherical]]]^2 - 
              probdensity, {x, 0.0001, 
              5}, {y, -0.0001, -5}, {z, -0.0001, -5}, 
            MaxRecursion -> 2, 
            DisplayFunction -> Identity]; 
Show[nwu, 
          nwd, swu, swd, neu, ned, seu, sed, 
          DisplayFunction -> $DisplayFunction]]

Find an estimate of the maximum probability density

pdens = 
      NMaximize[
        Abs[HAtom[2, 1, 1, 
              CoordinatesFromCartesian[{x, y, z}, Spherical]]]^2, {x, y, 
          z}]

It takes about a minute to compute and display one probability density contour

H221[0.9]  =  HAtomShow[2,  1, 1, 0.9*pdens[[1]]]

pdens = 
      NMaximize[
        Abs[HAtom[3, 2, 1, 
              CoordinatesFromCartesian[{x, y, z}, Spherical]]]^2, {x, y, 
          z}]

The next one takes a couple minutes, presumably because of the small probability values?

H321[0.5]  =  HAtomShow[3, 2, 1, 0.5*pdens[[1]]]

Show[H321[0.5], 
      ViewPoint -> {1.978,  -0.028,  -0.294}]