# HW 8
# Ha LUU
# 2/21/2014

############################################################################
#Problem 1

#ApplyNCg(f,x,n,a,b): inputs an expression f in the variable x,
#and a positive integer n, 2 variable a,b as interval, applies the Newton-Cotes formula 
#found by FindNC(n)  to approximate the integral of f(x) from a to b
ApplyNCg:=proc(f,x,n,a,b) local L,i:
L:=FindNC(n):
(b-a)*(add(L[i]*subs(x=(i-1)*(b-a)/n+a,f),i=1..n+1)):
end:

##############################
#Problem 2

ApplyNCgg:= proc(f,x,n,a,b,K) local L,i,m,j:

for j from 0 to K do
	m[j]:=a+(j)*(b-a)/K :
	#print(m[j]):
od:

for i from 0 to K-1 do
	L[i]:=(ApplyNCg(f,x,n,m[i],m[i+1])):
	#print(L[i]):
od:

add(L[i],i=0..K-1):

end:
#################
#Problem 3
###################################################
#BestNC(1/(1+x^2), x, 32) : is the the pair [2,16]
#BestNC(1/(1+x^4),x, 32): is the pair [8,4]
#BestNC(exp(x), x,32) : is the pair [4,8]
#BestNC(log(1+x),x, 32) : is the pair[4,8]
###################################################

## Findfactor

Findfactor:=proc(M) local n,i,j,k,L:

j:=1:

for i from 1 to M do
	if type(M/i,`integer`)=true then
		L[j][1]:=i:
		L[j][2]:=M/i:
		j:=j+1:
	fi:
od:

k:=j-1:

[seq([seq(L[j][k],k=1..2)],j=1..k)]:

end: 

## BestNC(f,x,M)
BestNC:=proc(f,x,M)  local L,K,i,j,errval:

L:= Findfactor(M):

errval:=1:

j=1:

for i from 1 to nops(L) do
	K[i]:= abs(evalf(ApplyNCgg(f,x,L[i][1],0,1,L[i][2]) - int(f,x=0..1))):
	if K[i]<errval then
		j:=i:
		errval:=K[i]:
	fi:
	
od:

print([L[j][1],L[j][2]]);

end:

########################################################################
#Problem 4

ApplyGaussian:=proc(f,x,n) local i,L :

L:=FindGauss(n):

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

end:
#######################################################################################################################################################################
$Problem 5

#for f:=1/(1+x^2); n:=1,2,3,4
#evalf(ApplyGaussian(f,x,n)-int(f,x=0..1));   0.146018365e-1  0.14870823e-2  -0.1311283e-3    0.481277978325867e-5
#evalf(ApplyNC(f,x,n)-int(f,x=0..1));        -0.353981635e-1 -0.20648302e-2  -0.7827789e-3    0.1312483e-3
# The ApplyGaussian is much more accurate function yet at n=5, it takes too much time to evaluate that I have to stop the program


######################################################################################################

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:


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


#ApplyNC(f,x,n,a,b): 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:

######################################################################