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