Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

Sign up
Here's how it works:
  1. Anybody can ask a question
  2. Anybody can answer
  3. The best answers are voted up and rise to the top

I was using Mathematica to get the series solution for Legendre equation. But when I get the recurrence relation and use RSolve:

RSolve[Subscript[a, n + 2] == ((n - l) (n + l + 1))/((1 + n) (2 + n)) Subscript[a, n], 
   Subscript[a, n], n] // FullSimplify // TraditionalForm

(where l and n are both non-negative)

Just for calculating the coefficient of the series, I only find the solution is strange:

{{Subscript[a, n]->Γ(1/2 (-l+n))/Γ(-(l/2)) (2^n (Subscript[c, 2] (-1)^n+Subscript[c, 1]) Γ(1/2 (l+n+1)))/(Γ(n+1) Γ(1/2 (l+1)))}}

I was confused that when l and n are even and n < l,the answer will be ComplexInfinity. But it is so counter-intuitive!

share|improve this question
What happens after applying FullSimplify[]? – J. M. May 11 '15 at 8:00

I think you need to use this form (subscripts removed because they're messy):

FullSimplify[RSolve[a[n + 2] == ((n + k + 1) (n + k) - 
 l (l + 1))/((n + k + 1) (n + k + 2)) a[n], a[n], n]]

(* you should get...*)
result = (2^n (C[1]+(-1)^n C[2])Gamma[1+k]Gamma[1/2 (k-l+n)]Gamma[
 1/2 (1+k+l+n)])/(Gamma[(k-l)/2] Gamma[1/2 (1+k+l)]Gamma[1+k+n])

Using this, you can now set k=0 to get the even half of the coefficients, and k=1 to get the odd half. This is the reason you were getting weird results, because you were trying to calculate the other half of coefficients from the wrong recurrence relation. NOTE: there are generated parameters C[1], C[2] because you haven't got specific initial conditions.

Now I think you can do:

evenCoeff = Table[FullSimplify[result /. {k -> 0}], {n, 0, 6, 2}]
oddCoeff = Table[FullSimplify[result /. {k -> 1}], {n, 0, 6, 2}]

Look at equations (4.13) in this text for more. I'm afraid I've never covered the Legendre equation or its polynomials so I'm only going on what I've read. The page on also says this which is suggestive of you get undefined results if you try to look for the even coefficients with the odd reccurence relation and vice versa:

If l is an even integer, the series y_1(x) reduces to a polynomial of degree l with only even powers of x and the series y_2(x) diverges. If l is an odd integer, the series y_2(x) reduces to a polynomial of degree l with only odd powers of x and the series y_1(x) diverges.

share|improve this answer
Yes,I just think that Gamma may appear oddly when n is not necessary positive. – WateSoyan May 11 '15 at 10:18
I'm still not sure about the correctness of my answer but I think that's approximately the reason why it was giving complex infinity. Gamma returns complex infinity on all non-positive integers as well. You could try FullSimplify[... , Element[n|l,Integers]&&n>0&&l>0] see if that helps. – Histograms May 11 '15 at 10:35

Note function LegendreP gives normalized orthogonal polynomials.

Approaching your recursion by nesting can reproduce these (though not proving where coefficients vanish). So for illustrative purposes (changing l to j):

f[n_, j_] := - (j + n + 1) (j - n)/((n + 1) (n + 2));
es[a_, j_, n_] := 
 First@Nest[{f[#[[2]], j] #[[1]], #[[2]] + 2} &, {a, 0}, n/2]
os[a_, j_, n_] := 
  First@Nest[{f[#[[2]], j] #[[1]], #[[2]] + 2} &, {a, 1}, Floor[n/2]];
esol[x_, j_?EvenQ] := (es[1, j, #] & /@ Range[0, 2 j, 2]).Table[
    x^(2 k), {k, 0, j}];
osol[x_, j_?OddQ] := (os[1, j, #] & /@ Range[1, 2 j + 1, 2]).Table[
    x^(2 k + 1), {k, 0, j}];
lp[n_?OddQ, x] := With[{c = osol[1, n]}, Expand[osol[x, n]/c]]
lp[n_?EvenQ, x] := With[{c = esol[1, n]}, Expand[esol[x, n]/c]]


Grid[{Expand@LegendreP[#, x], lp[#, x]} & /@ Range[8], Frame -> All]

enter image description here

share|improve this answer

Your Answer


By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.