# legendre.mpl -- Compute Legendre polynomials and their roots
#
# This Maple script creates the file legendre.m containing the
# roots and weights for Gaussing quadrature in the format used 
# in the function gauss.m
#
# Computations are done in quad precision and then rounded to
# double precision to ensure accuracy.
#
# Written Fall 2005 by Eric Olson for Math/CS 466/666

restart;
kernelopts(printbytes=false):
with(LinearAlgebra):
Digits:=32:
n:=4:
for k from 0 to n
do
  p[k]:=x^k;
  for j from 0 to k-1
  do
    p[k]:=p[k]-p[j]*int(p[k]*p[j],x=-1..1);
  end:
  p[k]:=expand(p[k]/sqrt(int(p[k]^2,x=-1..1)));
end:
for k from 1 to n
do
  r[k]:=Vector([fsolve(p[k]=0,x=-1..1)]);
end:
for k from 1 to n
do
  V:=VandermondeMatrix(r[k]);
  b:=Vector([seq(int(x^k,x=-1..1),k=0..k-1)]);
  w[k]:=LinearSolve(Transpose(V),b);
end:
for k from 1 to n
do
  G[k]:=unapply(DotProduct(w[k],Map(f,Vector(r[k]))),f);
end:
with(CodeGeneration):
# Warning, the protected name Matlab has been redefined and unprotected
rn:=Vector([seq(r[k],k=1..n)]):
wn:=Vector([seq(w[k],k=1..n)]):
fd:=fopen("legendre.m",WRITE): 
fprintf(fd,"roots=["):
for k from 1 to Dimension(rn)
do
  fprintf(fd," %.16g",rn[k]);
end:
fprintf(fd,"];\nweights=["):
for k from 1 to Dimension(wn)
do
  fprintf(fd," %.16g",wn[k]);
end:
fprintf(fd,"];\n"):
close(fd):