# Constructing recurrence equation for Jacobi series coefficients # of a solution of a linear differential equation. # ################################################################### # # Version 1.0 for MapleV release 4 # # Pawel Wozny (Pawel.Wozny@ii.uni.wroc.pl), 07/02/2001 # ################################################################### # # Let f(x) be a solution of a linear differential equation # # n # ----- # \ i # (1) ) v_i(x) D f(x) = g(x), # / # ----- # i = 0 # d # where v_i(x) are polynomials, and D:= --. # dx # Let a_k[f] be the Fourier coefficients of the function f(x) in the # expansion # ----- # \ # (2) f(x) = ) a_k[f] P_k(x), # / # ----- # k # # where P_k(x) is the k-th monic 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 # ~~~~~~~~~~~~~~ # - Jacobi(V) # # Input # ~~~~~ # - vector V := [v_0(x),v_1(x),...,v_n(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 (see (*)) is given on the screen such that # # B(k)= w*T(a_k[g]) # ################################################################### # # Example: # ~~~~~~~~ # Let f(x)=x*exp(x). This function satisfies the first-order differential # equation # # xf'(x)-(x+1)f(x)=0. # # We are loking for the recurrence equation for a_k[f]. We should initialize # Maple session as follows: # # >read(`jacobi.mpl`): # >Jacobi([-x-1,x]); # ################################################################### with(LREtools): bD:=f->D(unapply(f,x))(x): bV:=f->(x-zeta)*bD(f)+ni(zeta)*f: bK:=f->bV(bD(f)): bU:=f->sigma(x)*bD(f)+tau(x)*f: bL:=f->bU(bD(f)): zeta:=1: zeta2:=-1: _gamma:=alpha+beta+1: _delta:=alpha-beta: ni:=z->if z=1 then alpha+1 else beta+1 fi: mi:=(q,r,z)->pochhammer(k+ni(z),q+r)/(product(delta0(k+v),v=1..q)): sigma:=x->x^2-1: tau:=x->(_gamma+1)*x+_delta: h:=k->2^(2*k+_gamma)*k!*GAMMA(k+alpha+1)*GAMMA(k+beta+1)/ (GAMMA(2*k+_gamma+1)*pochhammer(k+_gamma,k)): lambda:=k->-k*(k+_gamma): ksi0:=k->4*k*(k+alpha)*(k+beta)*(k+_gamma-1)/(pochhammer(2*k+_gamma-2,3)* (2*k+_gamma-1)): ksi1:=k->-_delta*(_gamma-1)/((2*k+_gamma-1)*(2*k+_gamma+1)): ksi2:=k->1: delta0:=k->4*(k+alpha)*(k+beta)*pochhammer(k+_gamma-1,2)/ (pochhammer(2*k+_gamma-2,3)*(k+_gamma)*(2*k+_gamma-1)): delta1:=k->2*_delta/((2*k+_gamma-1)*(2*k+_gamma+1)): delta2:=k->-1/(k+_gamma): pX:=f->ksi0(k)*shift(f,k,-1)+ksi1(k)*f+ksi2(k)*shift(f,k,1): pD:=f->delta0(k)*shift(f,k,-1)+delta1(k)*f+delta2(k)*shift(f,k,1): lam_pI:=f->-lambda(k)*f: lambda_pD:=f->-lambda(k)*pD(f): pP:=f->f+zeta/2*pochhammer(2*k+_gamma,2)*pochhammer(2*k+_gamma+1,2)/ ((k+ni(zeta2))*(k+_gamma)*pochhammer(2*k+_gamma+m,2))*shift(f,k,1): pR:=f->delta0(k)*shift(f,k,-1)-2*zeta*(k+ni(zeta2)+m-1)/ pochhammer(2*k+_gamma+m-1,2)*f: pS:=proc(z,i,j,f) global m,zeta,zeta2; local l,S; S:=(pP@@0)(f); if z=1 then zeta:=z; zeta2:=-1 else zeta:=z; zeta2:=1 fi; for l from j to i do m:=l; S:=pP(S) od; S end: pZ:=proc(z,n,f) pS(z,n,1,f) end: pM:=proc(z,i,j,f) global m,zeta,zeta2; local l,R; R:=(pR@@0)(f); if z=1 then zeta:=z; zeta2:=-1 else zeta:=z; zeta2:=1 fi; for l from j to i by -1 do m:=l; R:=pR(R) od; R end: pN:=proc(z,n,f) pM(z,1,n,f) end: find_bQi:=proc(i,w) local wx,a,b,d,p,r,s,m,omega,eps; if simplify(w)=0 then RETURN([0,0,0,0,0,0,0]) fi; wx:=unapply(simplify(w),x); d:=degree(wx(x),x); a:=degree(gcd((x-1)^d,wx(x)),x); b:=degree(gcd((x+1)^d,wx(x)),x); p:=simplify(wx(x)/((x-1)^a*(x+1)^b)); r:=min(a,b); s:=max(a,b)-r; m:=floor((i+1)/2); omega:=i mod 2; if r=a then eps:=-1 else eps:=1 fi; if r>=m then [0,0,eps,0,m-omega,omega,simplify(p*(x^2-1)^(r-m)*(x-eps)^s)] elif (r+s>=m and m>r) then [0,m-r-omega,eps,omega,r,0,simplify(p*(x-eps)^(r+s-m))] else [i-2*r-2*s,s,eps,0,r,0,p] fi end: canonical_form:=proc(n,p) global zeta; local i,j,k,w,bQ,q,P; P[n]:=p; for k from 0 to n do w[k]:=coeff(P[n],(bD@@k)(f(x)),1) od; bQ:=array(0..n); for i from n to 1 by -1 do q:=find_bQi(i,w[i]); bQ[i]:=q; zeta:=q[3]; P[i-1]:=P[i]-(bD@@q[1]@bK@@q[2]@bV@@q[4]@bL@@q[5]@bU@@q[6])(q[7]*f(x)); for j from 0 to i-1 do w[j]:=coeff(P[i-1],(bD@@j)(f(x)),1) od od; bQ[0]:=[0,0,0,0,0,0,factor(simplify(w[0]))]; bQ end: poly_to_oper:=proc(w,f) local j,i,p; 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: integral_canonical_form:=proc(n,Q) local Omega,z,z2,t,i,j,l,p,p2,q,q2,e,epsilon,epsilon2, pom,O,sO,oO,stp,stm,sm,sp; z:=1; z2:=-1; Omega[z]:=[]; Omega[z2]:=[]; stp:=0; stm:=0; p[z]:=0; p[z2]:=0; q[z]:=0; q[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]); q[z]:=max(q[z],Q[i][1]+2*Q[i][2]+Q[i][4]); else Omega[z2]:=[op(Omega[z2]),i]; p[z2]:=max(p[z2],Q[i][1]); q[z2]:=max(q[z2],Q[i][1]+2*Q[i][2]+Q[i][4]) fi od; q[z]:=q[z]-p[z]; q[z2]:=q[z2]-p[z2]; if is(p[z]+q[z]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)); o1:=oper[4][1]; o2:=oper[4][2]; o3:=oper[4][3]; 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; print('w'=shift(w,k,prze)); print('s'=prze,'omega'=o1,'d'=o2,'e'=o3); print(T=1/'h'(k+prze)*E^'s'*Z['omega']^`(d)`*D^'e'); [sm+prze,sp+prze,A,shift(w,k,prze)] end: