Mathematica to Solve Legendre Differential Equation

<< DiscreteMath`RSolve`

First, define the infinite series that we assume the solution looks like (y = (∑_(i = 0)^∞) (Subscript[a, n]) 
            (x^n)). The Gf stands for generating function, and allows us to work with the infinite series symbolically.

Clear[α, c, a]

s[1] = PowerSum[a[n], {x, n, 0}]

Next, define the Legnedre differential equation:

L[y_] := (1 - x^2) D[y, {x, 2}] - 
        2  x D[y, x] + α (α + 1) y

Substitute our series solution into the differential equation:

s[2]  :=  L[s[1]]

Pull out the coefficient multiplying (x^n): (this does all the tedious series manipulations for us!)

s[3]  =  SeriesTerm[s[2], {x, 0, n}]

Solve the recurrence relation:

expr[1] = Solve[s[3] == 0, a[n + 2]]

From this equation, we can write down the first few terms and try to recognize a pattern. Here, we will get the first few terms in the two series solutions. From these, the pattern is easily recognized, although it can be inconvenient to write down. We will go out to a power num of x.

Clear[c, a]
num  =  9;
c[α, 0] = c1;
c[α, 1] = c2;
expr[2] = 
        Simplify[
          Part[expr[1], 1, 1, 2]  /.  
            n -> n - 2];
c[α_, n_]  :=  c[α, n] = 
        Simplify[
          expr[2] /. 
            a[n - 2] -> c[α, n - 2]]
TableForm[
        Table[{n, c[α, n]}, {n, 2, 
            num}]];
sol = Sum[c[α, n] x^n, {n, 0, num}]

Now, we can collect the two solutions:

sol  = Collect[sol, {c1, c2}];
sol1[x_, α_]  =  
      Coefficient[sol, c1]
sol2[x_, α_]  =  Coefficient[sol, c2]

From the above we search for a pattern. If we can't find a pattern, we can still generate the first few Legendre polynomials. Here is the pattern:

\(\(y\_1\) \((x)\) = \(1 - \(\(α \((1 + α)\)\)\/\(2!\)\) 
                x\^2 + \(\(α\ \((α - 
                          2)\) \((1 + α)\)\ \((3 + α)\)\)\/\(4!\)\) x\^4 -  ... \) + \(\((\(-1\))\)\^m\) \(\(\(\(α \((α - 
                              2)\) \((α - 4)\) ... \) \((α - 
                          2  m + 
                          2)\) \((1 + α)\) \((3 + α)\) ... \) \((α + 2  m - 1)\)\)\/\(\((2  m)\)!\)\) 
            x\^\(2  m\) -  ... 
\)
\(\(y\_2\) \((x)\) = \(x - \(\(\((α - 
                          1)\) \((2 + α)\)\)\/\(3!\)\) 
                x\^3 + \(\(\((α - 1)\) \((α - 
                          3)\) \((2 + α)\) \((4 + α)\)\)\/\(5!\)\) 
                x\^5 -  ... \) + \(\((\(-1\))\)\^m\) \(\(\(\(\((α - 
                              1)\) \((α - 3)\) ... \) \((α - 
                          2  m + 
                          1)\) \((2 + α)\) \((4 + α)\) ... \) \((α + 2  m)\)\)\/\(\((2  m + 1)\)!\)\) 
            x\^\(2  m + 1\) +  ... \)

Notice that for α zero or a positive even integer, the first series solution will eventually truncate, and for α odd the second series solution will eventually truncate. These truncated series form the Legendre polynomials.

LL[n_, x_]  :=  If[OddQ[n], sol2[x, n] , sol1[x, n] ]
TableForm[Table[{n, LL[n, x], LL[n, 1]}, {n, 0, num}]]

These still need to be standardized to have (Subscript[P, n]) (1) = 1to produce the Legendre polynomials you see in the literature. Compare with the actual Legendre polynomials, using Mathematica's intrinsic definition of the Legendre polynomial, below.

TableForm[
      Table[{n, LegendreP[n, x], LegendreP[n, 1]}, {n, 0, num}]]

Normalizations:

TableForm[
      Table[(∫_(-1)^1) LegendreP[n, x] 
            LegendreP[m, x] \[DifferentialD]x, {n, 0, num}, {m, 0, 
          num}]]

Here is a plot of the first few Legendre polynomials:

len = 11;
\(\(\(Do[
      plot[n] = 
        Plot[{LegendreP[n, x]}, {x, \(-1\), 1}, 
          Background -> GrayLevel[ .026], 
          PlotRange -> {{\(-1\), 1}, {\(-1\), 1}}, 
          PlotStyle -> {Hue[n/len]}, 
          DisplayFunction -> Identity], {n, 0, len}]\)\(
\)
    \)\)
Table[{plot[n], plot[n + 1], plot[n + 2]}, {n, 0, len, 
          3}];
Show[GraphicsArray[%], DisplayFunction -> $DisplayFunction, 
      ImageSize -> {650, Automatic}]