(* legendre.m math < legendre.m Legendre polynomials P_n(x) = ( p[[n+1]] in the program ) 2005.11.15 *) Clear[x]; main[]:=Module[{}, nn=10; (* Definition *) p={}; For[n=0,n<=nn,n++, f=Expand[D[(x^2-1)^n,{x,n}]/((2^n)*(n!))]; p=Append[p,f]; ]; (* the recursion formula *) q={1,x}; For[n=2,n<=nn,n++, g=Expand[((2*n-1)*x*q[[n]]-(n-1)*q[[n-1]])/n]; q=Append[q,g]; ]; (* check *) For[n=0,n<=nn,n++, If[!SameQ[p[[n+1]],q[[n+1]]], Print["ERROR !!"]; Print["Definition => P",n," = ",InputForm[p[[n+1]]]]; Print["Recursion Formula => P",n," = ",InputForm[q[[n+1]]]]; Exit[]; ]; ]; (* roots *) rlist={}; For[n=0,n<=nn,n++, rlist=Append[rlist,(x /. NSolve[p[[n+1]]==0,x])]; ]; (* Coeeicients c_n *) clist={}; For[n=0,n<=nn,n++, t=rlist[[n+1]]; c={}; For[k=1,k<=n,k++, h=1; For[j=1,j<=n,j++, If[j!=k, h*=(x-t[[j]])/(t[[k]]-t[[j]])]; ]; c=Append[c,Integrate[h,{x,-1,1}]]; ]; clist=Append[clist,c]; ]; (* Output *) Print[]; Print[]; For[n=0,n<=nn,n++, Print[]; Print["P",n," = ",InputForm[p[[n+1]]]]; If[n>0, Print["{t_k} : ",InputForm[rlist[[n+1]]]]; Print["{c_k} : ",InputForm[clist[[n+1]]]]; ]; ]; (* inner product For[m=0,m<=nn,m++, For[n=m,n<=nn,n++, ip=Integrate[p[[m+1]]*p[[n+1]],{x,-1,1}]; Print["< P",m," , P",n," > = ",InputForm[ip]]; ]; ]; *) ]; main[];