#C9.txt: Feb. 20, 2014, motivating orthogonal polynomials
#and rediscovering Gaussian quadrature

HelpOld:=proc(): print(`InterPolF(f,x,L) , FindGauss(n) `): 
print(`GuessRF(L,n) `):
end:

Help:=proc(): print(`ApI(f,x,X), Pn(n,x), FindGaussSmart(n) `): 
print(` ApplyQ(f,x,Q) `):
end:

#Old stuff needed today

#LagP(L,x): inputs a list L=[x1,...,xn]
# of numbers (or symbols)
#and a variable name x, and outputs
#the poly (x-x1)*(x-x2)* ...*x(x-xn)
LagP:=proc(L,x) local i:
mul(x-L[i],i=1..nops(L)):
end:
#LIF(L,x): inputs a list of pairs L
#L=[[x1,y1], ..., [xn,yn]] and outputs
#the unique polynomial in x of degree <n
#such that P(x1)=y1, ..., P(xn)=yn
LIF:=proc(L,x) local i,Ross,s, Sowmya:

Ross:=LagP([seq(L[i][1],i=1..nops(L))],x):

s:=0:

for i from 1 to nops(L) do
Sowmya:=normal(Ross/(x-L[i][1])):
#Sowmya:=Sowmya/subs(x=L[i][1],Sowmya):
Sowmya:=Sowmya/subs(x=L[i][1],diff(Ross,x)):
s:=s+Sowmya*L[i][2]:

od:

s:
end:

#InterPol(L,x):
#inputs a list of pairs [[x1,y1],[x2,y2], ..., [xn,yn]] of numbers (or symbols)! 
#and a variable name x,and outputs the unique polynomial of degree <= n-1 , 
#let's call it P(x), such that
#P(x1)=y1, P(x2)=y2, ..., P(xn)=yn,
#For example, InterPol([1,3],x); should give 3.
InterPol:=proc(L,x) local P,n,i,a,var,eq:

if not (type(L,list) and {seq(type(L[i],list),i=1..nops(L))}={true} and {seq(nops(L[i])-2,i=1..nops(L))}={0} ) then
print(`The first argument`, L, `should be a list of pairs`):
RETURN(FAIL):
fi:

n:=nops(L):

P:=add(a[i]*x^i,i=0..n-1):
var:={seq(a[i],i=0..n-1)}:

eq:={seq(numer(subs(x=L[i][1],P)-L[i][2])=0,i=1..nops(L))}:

var:=solve(eq,var):

if var=NULL then
 RETURN(FAIL):
else
 RETURN(subs(var,P)):
fi:

end:
#InterPolF(f,x,L):
#inputs an expression f, in the variable x, representing
#a function, and a list L=[x1,...,xn], of numbers (or symbols)
#outputs the unique polynomial of degree<= n-1, 
#let's call it P(x), such that
#P(x1)=f(x1), ..., P(xn)=f(xn). For example,
#InterPolF(x^3,[1,2,3,4],x); should return x^3
InterPolF:=proc(f,x,L) local i:
InterPol([seq([L[i],subs(x=L[i],f)], i=1..nops(L))], x):
end:

#FindGauss(n): inputs a pos. integer n
#and finds, using non-linear algebra
#the Gaussian Quadrature numerical (appx.)
#integration with clever
#intermetiate points [x1,..., xn]
#for int(f(x),x=0...1) that is EXACT
#for polynomials of degree <=2*n-1
#int(f(x),x=0..1)= H[1]*f(x1)+ ...
#+ H[n]*f(xn). The output is the pair
#[[x1,.., xn], [H[1], ..., H[n]]:
FindGauss:=proc(n) local x,H,var,eq,i:
option remember:


var:={seq(H[i],i=1..n), seq(x[i],i=1..n)}:

eq:= {seq( add(H[i]*x[i]^j,i=1..n)-1/(j+1), j=0..2*n-1)}:

var:=solve(eq,var):

[subs(var,[seq(x[i],i=1..n)]),
subs( var, [seq(H[i],i=1..n)])]:

end:



#GuessRF1(L,d,n): inputs a list of numbers (or expressions)
#L, a non-neg. integer d, and a (discrete) variable name
#n and tries to guess a rational function of n of degree
#d, such that L[i]-R(i)=0 for all i from 1 to nops(L)
GuessRF1:=proc(L,d,n) local R,a,b,i,j,eq,var:

if nops(L)<=2*d+6 then
 print(`Make list bigger`):
 RETURN(FAIL):
fi:

R:=(n^d+add(a[i]*n^i,i=0..d-1))/add(b[j]*n^j,j=0..d):

var:={seq(a[i],i=0..d-1), seq(b[j],j=0..d)}:

eq:={seq( numer(L[i]-subs(n=i,R)) =0,i=1..2*d+6)}:

var:=solve(eq,var):

if var=NULL then
 RETURN(FAIL):
fi:

R:=subs(var , R):


if {seq( numer(L[i]-subs(n=i,R)),i=1..nops(L))}<>{0}  then
 print(`It worked up to`, 2*d+6, `terms, but that's life `):
 FAIL:
else
 R:
fi:


end:


#GuessRF(L,n): inputs a list of numbers (or expressions)
#L, and a (discrete) variable name
#n and tries to guess a rational function of n of degree
#<=(nops(L)-6)/2 such that L[i]-R(i)=0 for all i from 1 to nops(L)
GuessRF:=proc(L,n) local d, katie:

for d from 0 to (nops(L)-6)/2 do
  katie:=GuessRF1(L,d,n):

   if katie<>FAIL then
      RETURN(katie):
   fi:

od:

FAIL:

end:


###END OF OLD stuff


#inputs an expression f of x (describing a function f of X)
#and finds an "approximation" of int(f,x=0..1);
#by first finding the Lagrange interpolating polynomial
#of f(x) at the list of points X=[x1,...] and then
#intergating the polynomial:
ApI:=proc(f,x,X) 
int(InterPolF(f,x,X),x=0..1):

end:

#Pn(n,x): the unique momic polynomial of x of degree n
#such that int(Pn(n,x)*x^i,x=0..1)=0 for i=0..n-1
Pn:=proc(n,x) local eq,var,i,P,a:
option remember:
P:=x^n+add(a[i]*x^i,i=0..n-1):

var:={seq(a[i],i=0..n-1)}:

eq:={ seq( int(P*x^i,x=0..1)=0, i=0..n-1)}:

sort(subs(solve(eq,var),P)):
end:

#FindGaussSmart(n): Smart version of FindGauss(n) 
#inputs a pos. integer n
#and finds, using non-linear algebra
#the Gaussian Quadrature numerical (appx.)
#integration with clever
#intermetiate points [x1,..., xn]
#for int(f(x),x=0...1) that is EXACT
#for polynomials of degree <=2*n-1
#int(f(x),x=0..1)= H[1]*f(x1)+ ...
#+ H[n]*f(xn). The output is the pair
#[[x1,.., xn], [H[1], ..., H[n]]:
FindGaussSmart:=proc(n) local X,x,H:
option remember:

X:= [fsolve(Pn(n,x))]:

H:=[seq(int(LIF([seq([X[i1],0],i1=1..i-1),[X[i],1],
seq([X[i1],0],i1=i+1..n)   ],x),x=0..1),i=1..n)]:

[X,H]:


end:

#ApplyQ(f,x,Q)
#inputs a function f of x, and a quadrature rule
#q=[X,H], computes it 
ApplyQ:=proc(f,x,Q) local i:

add(Q[2][i]*subs(x=Q[1][i],f), i=1..nops(Q[1])):


end: