#C8.txt: Feb. 17, 2014. Gaussian quadrature and

Help:=proc(): print(`LagP(L,x) , LIF(L,x), FindNC(n) `):
print(`ApplyNC(f,x,n) , FindGauss(n) `):
end:

###stuff from C7.txt
#C7.txt (virtual class, distance learning)
Help7:=proc(): print(`InterPol(L,x), InterPolF(f,x,L) `): 
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:

##end of C7.txt


#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:


#FindNC(n): inputs a pos. integer n
#and finds, using linear algebra
#the Newton-Cotes numerical (appx.)
#intermetiate points [x1,..., xn]

FindNC:=proc(n) local H,var,eq,i:
option remember:

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

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

var:=solve(eq,var):

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

end:


#ApplyNC(f,x,n): inputs an expression f in the variable x,
#and a positive integer n, applies the Newton-Cotes formula 
#found by FindNC(n)  to approximate the integral of f(x) from 0 to 1
ApplyNC:=proc(f,x,n) local L,i:
L:=FindNC(n):
add(L[i]*subs(x=(i-1)/n,f),i=1..n+1):
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: