# Constructing recurrence equation for Meixner series coefficients # of a solution of a linear difference equation. # ################################################################### # # Version 1.01 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 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 # # (D+)(f(x)):=f(x+1)-f(x), (D-)(f(x)):=f(x)-f(x-1). # # Let a_k[f] be the Fourier coefficients of the function f(x) in the # expansion # ----- # \ # (2) f(x) = ) a_k[f] M_k(x), # / # ----- # k # # where M_k(x) is the k-th monic Meixner 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 # ~~~~~~~~~~~~~~ # - Meixner(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,w], where A[-r..r+d] is a table of A_l(k) # (l=-r..r+d) # - operator T such that B(k)= w*T(a_k[g]) # (see [1] Algorithm 3.1, Step 6, Eq. (**)) is given on the screen # ################################################################### # # Example: # ~~~~~~~~ # Let f(x)=product(x+i,i=0..n-1) (Pochhammer's polynomial). # This function satisfies the first-order # difference equation # # (x+n-1)(D+)(f(x))-nf(x)=0. # # We are loking for the recurrence equation for a_k[f]. We should initialize # Maple session as follows: # # >read(`meixner.mpl`): # >Meixner([],[-n,x+n-1]); # ################################################################### # # References # ~~~~~~~~~~ # [1] S. Lewanowicz, P. Wozny, Algorithms for construction of # recurrence relations for the coefficients of expansions in # series of classical orthogonal polynomials, Institute of # Computer Science, University of Wroclaw, Wroclaw 2001. # # http://www.ii.uni.wroc.pl/~pwo/programs.html # http://www.ii.uni.wroc.pl/~sle/publ/ # ################################################################### with(LREtools): k:='k': m:='m': bD:=f->if R=`+` then delta(f,x) else f-shift(f,x,-1) fi: bU:=f->sigma(x)*bD(f,x)+tau(x)*f: zeta_p:=-beta: zeta_m:=0: sigma:=x->if R=`+` then c*(beta+x) else x fi: tau:=x->beta*c+(c-1)*x: h:=k->k!*pochhammer(beta,k)*c^k/(1-c)^(beta+2*k): lambda:=k->(1-c)*k: ksi0:=k->c*k*(k+beta-1)/(1-c)^2: ksi1:=k->((c+1)*k+beta*c)/(1-c): ksi2:=k->1: delta_p0:=k->c*(1-beta-k)/((c-1)*(1-c)): delta_p1:=k->c/(1-c): delta_p2:=k->0: delta_m0:=k->delta_p0(k): delta_m1:=k->1+delta_p1(k): delta_m2:=k->delta_p2(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)*(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: szukaj_bQi:=proc(i,w,r) global R,zeta_p,zeta_m; local test1,j,wx,q,a,z1,z2,s; R:=r; if simplify(w)=0 then RETURN([0,0,0,r]) fi; test1:=true; a:=0; wx:=unapply(simplify(w),x); if r=`+` then z1:=zeta_p; s:=-1 else z1:=zeta_m; s:=1 fi; for j from 1 to i while test1 do if simplify(wx(-s*j+z1+s*i))=0 then a:=a+1 else test1:=false fi od; q:=simplify(wx(x)/product(sigma(x-s*i+s*k),k=1..a)); [i-a,a,(-s)^i*factor(simplify(shift(q,x,s*i))),r] end: postac_kanoniczna:=proc(n,p,r) global R; local i,j,k,w,bQ,q,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(collect(P[n],f(x+s*k)),f(x+s*k),1) od; bQ:=array(0..n); for i from n to 1 by -1 do q:=szukaj_bQi(i,w[i],r); bQ[i]:=q; P[i-1]:=P[i]-(bD@@q[1]@bU@@q[2])(q[3]*f(x)); for j from 0 to i-1 do w[j]:=coeff(collect(P[i-1],f(x+s*j)),f(x+s*j),1) od od; bQ[0]:=[0,0,factor(simplify(w[0])),r]; bQ end: w_na_op:=proc(w,f) local i,p,j,ww; ww:=expand(w); for i from 0 to degree(ww,x) do p[i]:=coeff(ww,x,i)*(pX@@i)(f) od; sum(p[j],j=0..degree(ww,x)) end: calkuj_postac_kanoniczna:=proc(n,Q,r) global R; local p,i,j,t,oO,O,sO,stm,stp,l; p:=0; stm:=0: stp:=0; R:=r; for i from 0 to n do p:=max(p,Q[i][1]) od; for j from 0 to n do O[j]:=(pD@@(p-Q[j][1])@lambda_pD@@Q[j][2])(w_na_op(Q[j][3],b[k])); stm:=max(stm,degree(Q[j][3],x)+Q[j][2]+p-Q[j][1]); stp:=max(stp,degree(Q[j][3],x)) od; oO:=sum(O[t],t=0..n); sO:=0; for l from -stm to stp do sO:=sO+factor(simplify(coeff(collect(oO,b[k+l]),b[k+l],1)))*b[k+l] od; [p,sO,stm,stp] end: p_i_m:=proc(pp,pMp,stmp,stpp,pm,pMm,stmm,stpm) global R; local pUp,pUm,stm,stp,oO,sO,l; R:=`-`; pUp:=(pD@@pm)(pMp); R:=`+`; pUm:=(pD@@pp)(pMm); stp:=max(stpp,stpm); stm:=max(stmp+pm,stmm+pp); oO:=pUp+pUm; sO:=0; for l from -stm to stp do sO:=sO+factor(simplify(coeff(collect(oO,b[k+l]),b[k+l],1)))*b[k+l] od; [sO,stm,stp,0] end: MeixnerE:=proc(lm,lp) options `Copyright 2000 by Pawel Wozny`; global R; local Pp,Pm,np,nm,Kp,Km,m,p,i,j,n,l,oper,d,d2, A,B,pom,sp,sm,prze,o,w,O,O1,O2; nm:=nops(lm); np:=nops(lp); R:=`-`; if nm=0 then Pm:=0 else nm:=nm-1; Pm:=sum(simplify(lm[n+1])*f(x-n),n=0..nm) fi; R:=`+`; if np=0 then Pp:=0 else np:=np-1; Pp:=sum(simplify(lp[l+1])*f(x+l),l=0..np) fi; Km:=postac_kanoniczna(nm,Pm,`-`); Kp:=postac_kanoniczna(np,Pp,`+`); m:=calkuj_postac_kanoniczna(nm,Km,`-`); p:=calkuj_postac_kanoniczna(np,Kp,`+`); oper:=p_i_m(p[1],p[2],p[3],p[4],m[1],m[2],m[3],m[4]); sm:=0; sp:=0; d:=1; d2:=1; for i from -oper[2] to oper[3] do pom:=coeff(collect(oper[1],b[k+i]),b[k+i],1); B[i]:=factor(simplify(pom*h(k+i)/h(k))); if B[i]<>0 then d:=lcm(d,denom(B[i])); d2:=gcd(d2,numer(B[i])); sp:=max(sp,i); sm:=min(sm,i) fi od; prze:=-sm-floor((sp-sm)/2); w:=factor(simplify(d/d2)); o:=oper[4]; for j from sm+prze to sp+prze do A[j]:=0; if B[j-prze]<>0 then A[j]:=factor(simplify(shift(B[j-prze]*w,k,prze))); print(`A[`.j.`]`=A[j]) fi od; if prze=0 then O:=1 else O:='E^s' fi: if o[1]=0 then O1:=1 else O1:='(D^`+`)^u' fi: if o[2]=0 then O2:=1 else O2:='(D^`-`)^v' fi: print('w'=shift(w,k,prze)); print(`if`(prze=0,NULL,'s'=prze),`if`(o[1]=0,NULL,'u'=o[1]), `if`(o[2]=0,NULL,'v'=o[2])); print('T'='w'*1/'h'(k+prze)*O*O1*O2); [sm+prze,sp+prze,A,shift(w,k,prze)] end: Meixner:=proc(lm,lp) options `Copyright 2000 by Pawel Wozny`; global R; local Pp,Pm,np,nm,Kp,Km,m,p,i,j,n,l,oper,d,d2, A,B,pom,sp,sm,prze,o,w,O,O1,O2; nm:=nops(lm); np:=nops(lp); R:=`-`; if nm=0 then Pm:=0 else nm:=nm-1; Pm:=sum(simplify(lm[n+1])*(bD@@n)(f(x)),n=0..nm) fi; R:=`+`; if np=0 then Pp:=0 else np:=np-1; Pp:=sum(simplify(lp[l+1])*(bD@@l)(f(x)),l=0..np) fi; Km:=postac_kanoniczna(nm,Pm,`-`); Kp:=postac_kanoniczna(np,Pp,`+`); m:=calkuj_postac_kanoniczna(nm,Km,`-`); p:=calkuj_postac_kanoniczna(np,Kp,`+`); oper:=p_i_m(p[1],p[2],p[3],p[4],m[1],m[2],m[3],m[4]); sm:=0; sp:=0; d:=1; d2:=1; for i from -oper[2] to oper[3] do pom:=coeff(collect(oper[1],b[k+i]),b[k+i],1); B[i]:=factor(simplify(pom*h(k+i)/h(k))); if B[i]<>0 then d:=lcm(d,denom(B[i])); d2:=gcd(d2,numer(B[i])); sp:=max(sp,i); sm:=min(sm,i) fi od; prze:=-sm-floor((sp-sm)/2); w:=factor(simplify(d/d2)); o:=oper[4]; for j from sm+prze to sp+prze do A[j]:=0; if B[j-prze]<>0 then A[j]:=factor(simplify(shift(B[j-prze]*w,k,prze))); print(`A[`.j.`]`=A[j]) fi od; if prze=0 then O:=1 else O:='E^s' fi: if o[1]=0 then O1:=1 else O1:='(D^`+`)^u' fi: if o[2]=0 then O2:=1 else O2:='(D^`-`)^v' fi: print('w'=shift(w,k,prze)); print(`if`(prze=0,NULL,'s'=prze),`if`(o[1]=0,NULL,'u'=o[1]), `if`(o[2]=0,NULL,'v'=o[2])); print('T'='w'*1/'h'(k+prze)*O*O1*O2); print(o); [sm+prze,sp+prze,A,shift(w,k,prze)] end: