(* 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[];