################################################################### # # Warning!!! To run this program you need file qsum.mpl, which is available # ~~~~~~~~~~ on http://www.imn.htwk-leipzig.de/~koepf. # ################################################################### # # Constructing recurrence equation for the little q-Jacobi series # coefficients of a solution of a linear q-difference equation. # ################################################################### # # Version 1.0 for MapleV release 4 # # Copyright (c) Pawel Wozny (Pawel.Wozny@ii.uni.wroc.pl), 07/02/2001 # ################################################################### # # Let f(x) be a solution of a linear q-difference equation # # n m # ----- ----- # \ i \ j # (1) ) u_i(x) (D-) f(x) + ) v_j(x) (D+) f(x) = g(x), # / / # ----- ----- # i = 0 j = 0 # # where u_i(x), v_j(x) are polynomials, and # # f(q x) - f(x) # (D+) = -------------, # x (q - 1) # # f(x/q) - f(x) # (D-) = ------------- (0 < q <1 ). # x (1/q - 1) # # Let a_k[f] be the Fourier coefficients of the function f(x), # # ----- # \ # (2) f(x) = ) a_k[f] p_k(x), # / # ----- # k # # where p_k(x) is the k-th monic little q-Jacobi polynomial. This program # finds recurrence relation for a_k[f], of the form # # r + d # ----- # \ # (3) ) A_i(k) a_{i+r}[f] = B(k) # / # ----- # i = -r # # under assumption that a_k[g] are known explicitly. Here r is a natural # number and d î {0,1}. # # Main procedure # ~~~~~~~~~~~~~~ # - l_qJacobi(U,V) # # Input # ~~~~~ # - vectors U := [u_0(x),u_1(x),...,u_n(x)] and # V := [v_0(x),v_1(x),...,v_m(x)] # # Output # ~~~~~~ # - vector [-r,r+d,A], where A[-r..r+d] is a table of A_l(k) # (l=-r..r+d) # - operator T such that B(k)= T(a_k[g]) # (see Algorithm 3.1, Step 6, Eq. (**)) is given on the screen # ################################################################### # # Example: # ~~~~~~~~ # Let f(x)=x^n. This function satisfies the first-order # q-difference equation # # x(D+)(f(x))-[(q^n-1)/(q-1)]f(x)=0. # # We are loking for the recurrence equation for a_k[f]. We should initialize # Maple session as follows: # # >read(`l_jacobi.mpl`): # >l_qJacobi([],[-qnum(n),x]); # ################################################################### read(`d:/user/sle/classalg/qsum.mpl`): ##!!!!!! 14.04.03 with(LREtools): bE:=(f,i)->subs({x=q^i*x,X=q^i*X},f): qpochhammer:=(x,q,kk)->product(1-x*q^mm,mm=0..kk-1): qnum:=kk->(q^kk-1)/(q-1): eps:=()->if R=`+` then 1 else q fi: ww:=(n,q)->if (is(n=0) or is(q=1)) then 1 else if R=`+` then 1/(q^((n+1)*(n+2)/2-2*n-1)*(q-1)^n) else 1/(q^(-(n+1)*(n+2)/2+2*n+1)*(q^(-1)-1)^n) fi fi: bD:=f->if R=`+` then (bE(f,1)-f)/(X*(q-1)) else (bE(f,-1)-f)/(X*(q^(-1)-1)) fi: bU:=f->sigma(x)/eps()*bD(f,x)+tau(x)*f: bV:=f->if R=`+` then (x-zeta_p)*bD(f,x) else (x-zeta_m)/q*bD(f) fi: zeta_p2:=0: zeta_p:=1/(b*q): zeta_m2:=1: zeta_m:=0: ni_p:=(zp2,zm2)->qsimplify(2*tau(zm2)/(2*a*b*q^2*(zm2-zp2))): ni_m:=(zp2,zm2)->qsimplify(2*tau(zp2)/(2*(zp2-zm2))): gamma_p:=(zp2,zm2,k)->qsimplify(-q^(k-1)*lambda(k)/RR(k-1)+ni_p(zp2,zm2)): gamma_m:=(zp2,zm2,k)->qsimplify(1/(q-1-q^(1-k)*RR(k-1)/lambda(k))+ni_m(zp2,zm2)): sigma:=x->if R=`+` then a*q*x*(b*q*x-1) else x*(x-1) fi: sigmam:=x->x*(x-1): tau:=x->(a*q*(b*q*x-1)+1-x)/(q-1): lambda:=k->-(q^k-1)*(a*b*q^(k+1)-1)/(q^(k-1)*(q-1)^2): v_p:=k->qsimplify(-q^(k-1)*lambda(k)/RR(k-1)): v_m:=k->qsimplify(1/(q-1-q^(1-k)*RR(k-1)/lambda(k))): t_1:=diff(tau(x),x): t_0:=tau(0): s_2:=simplify(diff(sigmam(x),x,x)/2+(q-1)*t_1): s_1:=simplify(subs(x=0,diff(sigmam(x),x))+(q-1)*tau(0)): s_0:=sigmam(0): RR:=k->s_2*qnum(k)+t_1: h:=n->if n=0 then 1 elif n>0 then product(ksi0(k+ii),ii=1..n) else product(1/ksi0(k+ii),ii=n+1..0) fi: ksi0:=k->a*q^(2*k)*(q^k-1)*(a*q^k-1)*(b*q^k-1)*(a*b*q^k-1) /((a*b*q^(2*k)-1)^2*(a*b*q^(2*k)-q)*(a*b*q^(2*k+1)-1)): ksi1:=k->q^k*(1+a+a*(1+a)*b*q^(2*k+1)-a*(1+b)*q^k*(1+q)) /((a*b*q^(2*k)-1)*(a*b*q^(2*k+2)-1)): ksi2:=k->1: delta_p0:=k->qsimplify(-q^(1-k)*RR(k-1)*ksi0(k)/lambda(k)): delta_p1:=k->qsimplify(-q*t_0/((q+1)*lambda(k))- q*ksi1(k)*((q-1)*qnum(-k)*RR(k-1)+t_1)/ ((q+1)*lambda(k))): delta_p2:=k->qsimplify(s_2*qnum(k)/lambda(k)): delta_m0:=k->qsimplify(delta_p0(k)+(q-1)*ksi0(k)): delta_m1:=k->qsimplify(delta_p1(k)+(q-1)*ksi1(k)): delta_m2:=k->qsimplify(delta_p2(k)+(q-1)*ksi2(k)): pX:=f->ksi0(k)*shift(f,k,-1)+ksi1(k)*f+ksi2(k)*shift(f,k,1): pD:=f->if R=`+` then delta_p0(k)*shift(f,k,-1)+delta_p1(k)*f+delta_p2(k)*shift(f,k,1): else delta_m0(k)*shift(f,k,-1)+delta_m1(k)*f+delta_m2(k)*shift(f,k,1): fi: lambda_pD:=f->if R=`+` then -lambda(k)/q*(delta_m0(k)*shift(f,k,-1)+ delta_m1(k)*f+delta_m2(k)*shift(f,k,1)) else -lambda(k)*(delta_p0(k)*shift(f,k,-1)+ delta_p1(k)*f+delta_p2(k)*shift(f,k,1)) fi: pDelta:=f->if R=`+` then if m>0 then (pP@pR)(f)+(q-1)*(pQ@pR)(f)+zeta_p*(q-1)*f else delta_m0(k)*shift(f,k,-1)+delta_m1(k)*f+ delta_m2(k)*shift(f,k,1) fi else if m>0 then (pP@pR)(f)-(q-1)*(pQ@pR)(f)-zeta_m*(q-1)*f else delta_p0(k)*shift(f,k,-1)+delta_p1(k)*f+ delta_p2(k)*shift(f,k,1) fi fi: pP:=f->if R=`+` then f+ b*q^m*qpochhammer(a*b*q^(2*k+1),q,3)*(a*b*q^(2*k+2)-1)/ ((b*q^(k+1)-b*q*zeta_p2)*(1-zeta_p*a*b*q^(k+2))* (a*b*q^(k+1)-1)*qpochhammer(a*b*q^(2*k+m+1),q,2))*shift(f,k,1) else f+ q^(m-1)*qpochhammer(a*b*q^(2*k+1),q,3)*(a*b*q^(2*k+2)-1) /(q^k*(a*q^(k+1)-zeta_m)* (zeta_m2*b*q^(k+1)-1)*(a*b*q^(k+1)-1)* qpochhammer(a*b*q^(2*k+m+1),q,2))*shift(f,k,1) fi: pQ:=f->if R=`+` then v_p(k)*f+ qpochhammer(a*b*q^(2*k+1),q,3)*(a*b*q^(2*k+2)-1)*(1-a*b*q^(k+m+1))/ (a*q^(k+1)*(q-1)*(b*q^(k+1)-b*q*zeta_p2)*(1-zeta_p*a*b*q^(k+2))* (a*b*q^(k+1)-1)*qpochhammer(a*b*q^(2*k+m+1),q,2))*shift(f,k,1) else v_m(k)*f+ q^(m-1)*qpochhammer(a*b*q^(2*k+1),q,3)*(a*b*q^(2*k+2)-1)* (1-a*b*q^(k+m+1))/(q^k*(q-1)*(a*q^(k+1)-zeta_m)* (zeta_m2*b*q^(k+1)-1)*(a*b*q^(k+1)-1)* qpochhammer(a*b*q^(2*k+m+1),q,2))*shift(f,k,1) fi: pR:=f->if R=`+` then delta_p0(k)*shift(f,k,-1)+ a*q^k*(q-1)*(b*q^(k+m)-b*q*zeta_p2)*(1-zeta_p*a*b*q^(k+m+1))/ qpochhammer(a*b*q^(2*k+m),q,2)*f else delta_m0(k)*shift(f,k,-1)+ q^k*(q-1)*(a*q^(k+m)-zeta_m)*(zeta_m2*b*q^(k+m)-1)/ qpochhammer(a*b*q^(2*k+m),q,2)*f fi: pF:=f->if R=`+` then ni_p(zeta_p2,zeta_m2)*delta_m0(k)/gamma_p(zeta_p2,zeta_m2,k-1)* shift(f,k,-1)+ b*q^(k+1)*(q-1)*(a*q-zeta_m)*(1-zeta_p)* (zeta_m2*q^(k+m-1)-zeta_p2)/qpochhammer(a*b*q^(2*k+m),q,2)*f else q*ni_m(zeta_p2,zeta_m2)*delta_p0(k)/gamma_m(zeta_p2,zeta_m2,k-1)* shift(f,k,-1) -a*q^(2*k+m)*(b*q-zeta_m)*(q-1)* (zeta_m2*q^(k+m-1)-zeta_p2+zeta_m*(1+zeta_p2))*(a*b*q^2*zeta_p-1)/ qpochhammer(a*b*q^(2*k+m),q,2)*f fi: find_bQi:=proc(i,w,r) global R; local test1,test2,j,wx,wx2,wx3,qq,stala,aa,bb,z1,z2,s,m; R:=r; if qsimplify(w)=0 then RETURN([0,0,0,0,0,r]) fi; test1:=true; test2:=true; aa:=0; bb:=0; wx:=unapply(qsimplify(w),x); wx2:=unapply(qsimplify(w),x); wx3:=unapply(qsimplify(w),x); if r=`+` then z1:=0; z2:=1/(b*q); s:=-1; m:=a*b*q^2/eps() else z1:=1; z2:=0; s:=1; m:=1/eps() fi; for j from 1 to i while (test1 or test2) do if (test1 and simplify(wx2(z1*q^(s*(i-j))))=0) then aa:=aa+1; wx2:=unapply(qsimplify(wx2(x)/(x*q^(-s*(i-j))-z1)),x) else test1:=false fi; if (test2 and simplify(wx3(z2*q^(s*(i-j))))=0) then bb:=bb+1; wx3:=unapply(qsimplify(wx3(x)/(x*q^(-s*(i-j))-z2)),x) else test2:=false fi od; qq:=simplify(wx(x)/(product(x*q^(s*(k-i))-z1,k=1..aa)* product(x*q^(s*(l-i))-z2,l=1..bb))); if min(aa,bb)=bb then stala:=ww(i-aa,q)*ww(aa-bb,q)*ww(bb,q)*q^(s*((i-aa)*(aa-bb)+bb*(i-bb))); [i-aa,aa-bb,z1,bb,factor(qsimplify(eps()^(aa-bb)/(stala*m^bb)* subs(X=x,bE(qq*X^i,s*i)))),r] else stala:=ww(i-bb,q)*ww(bb-aa,q)*ww(aa,q)*q^(s*((i-bb)*(bb-aa)+aa*(i-aa))); [i-bb,bb-aa,z2,aa,factor(qsimplify(eps()^(bb-aa)/(stala*m^aa)* subs(X=x,bE(qq*X^i,s*i)))),r] fi end: canonical_form:=proc(n,p,r) global zeta_p,zeta_m,R; local i,j,k,w,bQ,qq,s,P; R:=r; P[n]:=p; if R=`+` then s:=1 else s:=-1 fi; for k from 0 to n do w[k]:=coeff(P[n],f(x*q^(s*k)),1) od; bQ:=array(0..n); for i from n to 1 by -1 do qq:=find_bQi(i,w[i],r); bQ[i]:=qq; if r=`+` then zeta_p:=qq[3] else zeta_m:=qq[3] fi; P[i-1]:=P[i]-(bD@@qq[1]@bV@@qq[2]@bU@@qq[4])(qq[5]*f(x)); for j from 0 to i-1 do w[j]:=coeff(P[i-1],f(x*q^(s*j)),1) od od; bQ[0]:=[0,0,0,0,factor(qsimplify(subs(X=x,w[0]))),r]; bQ end: poly_to_oper:=proc(w,f) local i,p,j; for i from 0 to degree(w,x) do p[i]:=coeff(w,x,i)*(pX@@i)(f) od; sum(p[j],j=0..degree(w,x)) end: pS:=proc(z,i,j,r,f) global R,m,zeta_p,zeta_p2,zeta_m,zeta_m2; local l,S; R:=r; S:=(pP@@0)(f); if r=`+` then if is(z=0) then zeta_p:=z; zeta_p2:=1/(b*q) else zeta_p:=z; zeta_p2:=0 fi else if is(z=1) then zeta_m:=z; zeta_m2:=0 else zeta_m:=z; zeta_m2:=1 fi fi; for l from j to i do m:=l; S:=pP(S) od; S end: pZ:=proc(z,n,r,f) pS(z,n,1,r,f) end: pU:=proc(z,n,r,f) global R,m,zeta_p,zeta_p2,zeta_m,zeta_m2; local l,Q; R:=r; Q:=(pQ@@0)(f); if r=`+` then if is(z=0) then zeta_p:=z; zeta_p2:=1/(b*q) else zeta_p:=z; zeta_p2:=0 fi else if is(z=1) then zeta_m:=z; zeta_m2:=0 else zeta_m:=z; zeta_m2:=1 fi fi; for l from 1 to n do m:=l; Q:=pQ(Q) od; Q end: pM:=proc(z,i,j,r,f) global R,m,zeta_p,zeta_p2,zeta_m,zeta_m2; local l,_R; R:=r; _R:=(pR@@0)(f); if r=`+` then if is(z=0) then zeta_p:=z; zeta_p2:=1/(b*q) else zeta_p:=z; zeta_p2:=0 fi else if is(z=1) then zeta_m:=z; zeta_m2:=0 else zeta_m:=z; zeta_m2:=1 fi fi; for l from j to i by -1 do m:=l; _R:=pR(_R) od; _R end: pN:=proc(z,n,r,f) pM(z,1,n,r,f) end: pW:=proc(zp,zm,i,j,r,f) global R,m,zeta_p,zeta_p2,zeta_m,zeta_m2; local l,V; R:=r; V:=(pF@@0)(f); if is(zp=0) then zeta_p:=zp; zeta_p2:=1/(b*q) else zeta_p:=zp; zeta_p2:=0 fi; if is(zm=0) then zeta_m:=zm; zeta_m2:=1 else zeta_m:=zm; zeta_m2:=0 fi; for l from j to i by -1 do m:=l; V:=pF(V) od; V end: pY:=proc(zp,zm,n,r,f) pW(zp,zm,1,n,r,f) end: integrate_canonical_form:=proc(n,Q,r) global R; local Omega,z,z2,i,j,l,p,p2,qq,q2,e,epsilon,epsilon2, pom,O,sO,oO,stp,stm,sm,sp,t; R:=r; if r=`+` then z:=0; z2:=1/(b*q) else z:=1; z2:=0 fi; Omega[z]:=[]; Omega[z2]:=[]; stp:=0; stm:=0; p[z]:=0; p[z2]:=0; qq[z]:=0; qq[z2]:=0; for i from 0 to n do if is(Q[i][3]=z) then Omega[z]:=[op(Omega[z]),i]; p[z]:=max(p[z],Q[i][1]); qq[z]:=max(qq[z],Q[i][1]+Q[i][2]); else Omega[z2]:=[op(Omega[z2]),i]; p[z2]:=max(p[z2],Q[i][1]); qq[z2]:=max(qq[z2],Q[i][1]+Q[i][2]) fi od; qq[z]:=qq[z]-p[z]; qq[z2]:=qq[z2]-p[z2]; if is(p[z]+qq[z]0 and j>0 and (zeta_p-zeta_m2)*(zeta_m-zeta_p2)=0) then if i0 then sp:=max(sp,i); sm:=min(sm,i) fi od; prze:=-sm-floor((sp-sm)/2); o:=oper[4]; for j from sm+prze to sp+prze do A[j]:=0; if B[j-prze]<>0 then A[j]:=factor(qsimplify(shift(B[j-prze],k,prze))); print(`A[`.j.`]`=A[j]) fi od; if o[1]=`+` then print('s'=prze,'zeta'=o[2],'eta'=o[6],'d'=o[3],'t'=o[4], 'u'=o[7],'v'=o[8]); print(T=1/'h'(k+prze)*E^'s'*Z['zeta']^`+(d)`*Y['zeta','eta']^`+(t)`* Z['zeta']^`+(t)`*(D^`-`)^'u'*(D^`+`)^'v'); else print('s'=prze,'epsilon'=o[2],'zeta'=o[5],'eta'=o[6],'d'=o[3], 't'=o[4],'u'=o[7],'v'=o[8]); print(T=1/'h'(k+prze)*E^'s'*Z['epsilon']^`-(d)`*Y['zeta','eta']^`+(t)`* Z['zeta']^`+(t)`*(D^`-`)^'u'*(D^`+`)^'v'); fi: [sm+prze,sp+prze,A] end: