#C10.txt, Feb. 24, 2014 Orthogonal polynomials in general Help:=proc(): print(`GuessRF(L,n) , AM(M,P,x)`): print(`PnM(M,n,x) , PnMseq(M,n,x), WtToM(w,n,x,a,b), AnM(M,n) `): end: #old stuff #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 old stuff #AM(M,P,x): inputs a sequence of moments x^i->M[i+1] #outputs the "integral" of polynomial P of x AM:=proc(M,P,x) local i: if nops(M)m PnM:=proc(M,n,x) local P,i,a,var,eq: option remember: P:=x^n+add(a[i]*x^i, i=0..n-1): var:={seq(a[i],i=0..n-1)}: #AM (M,x^i*P,x)=0, i=0, ...n-1 eq:={ seq(AM(M,x^i*P,x)=0, i=0..n-1)}: var:=solve(eq,var): if var=NULL then FAIL: else subs(var,P): fi: end: PnMseq:=proc(M,n,x) local i: [seq(PnM(M,i,x),i=0..n)]: end: #WtToM(w,n,x,a,b): the moments of f(x)->int(w(x)*f(x),x=a..b); WtToM:=proc(w,n,x,a,b) local i: [seq(int(w*x^i,x=a..b), i=0..n)]: end: #A(n)=Int(x*Pn(x)^2)/Int(Pn(x)^2) #B(n)=Int(x*Pn(x)*P_(n-1)(x))/Int(P_(n-1)^2) #AnM(M,n): the first n terms of the A(n) AnM:=proc(M,n) local i,andrew,x: andrew:=PnMseq(M,n,x): [seq(AM(M,x*andrew[i+1]^2,x)/AM(M,andrew[i+1]^2,x),i=1..n)]: end: