# $Source: /u/maple/research/lib/DEtools/src/RCS/RiemannPsols,v $ # $Notify: mvanhoei@daisy.uwaterloo.ca $ # # Solve the Riemann-Papperitz equation which is a Fuchsian equation # of order 2 with three (regular) singular points. # # Method: table lookup for the exponents, checking for integer differences, # and substituting the exponents in the solution. macro( e_int =`DEtools/kovacicsols/e_int`, ispolylinearODE =`dsolve/diffeq/ispolylinearODE`, sols =`DEtools/RiemannP_s`, sols2 =`DEtools/RiemannP_s2`, sols3 =`DEtools/RiemannP_s3`, S_and_E =`DEtools/RiemannP_SE`, sing =`DEtools/RiemannP_si` ): `DEtools/RiemannPsols` := proc(A,x) local ans,gg,AA,C,S,V; option `Copyright (c) 1997 Waterloo Maple Inc. All rights reserved. Author: M. van Hoeij`; if not type(A,list) then ans := DEtools[convertAlg](A,x); if not ans = FAIL and ispolylinearODE(ans[1],ans[2],op(x),'AA','gg')then if gg = 0 then return procname(AA,op(x)); else return []; fi; else return []; fi; fi; if not (_Env_odsolve_int::symbol and member(_Env_odsolve_int,['int,Int,`PDEtools/int`'])) then _Env_odsolve_int:=`PDEtools/int` fi: if nops(A)<>3 then [] elif indets([args],{radical,nonreal})<>{} then V := radfield(indets([args],{radical,nonreal})); eval(subs(V[2],procname(op(eval(subs(V[1],[args])))))) elif A[3]<>1 then # make it monic: procname([evala(Normal(A[1]/A[3])),evala(Normal(A[2]/A[3])),1],x,args) else # do some preprocessing: S:=-A[2]/2; # D -> D+S makes the coeff(..,D,1) vanish, and # gives the following coeff(..,D,0): C:=evala(A[1]-S^2+diff(S,x)); # Compute the finite singularities: V:=sing(C,args[2..nargs]); if V=FAIL then [] else sols([A[3],A[2],A[1]],S,C,V,x) fi fi end: sols:=proc(A,S,C,V,x) local i,D,a0,a1,d,s,v; option `Copyright (c) 1997 Waterloo Maple Inc. All rights reserved.`; D:=convert(V,`*`); d:=degree(D,x); # number of finite singularities a1:=collect(evala(Normal(A[2]*D)),x,evala); a0:=collect(evala(Normal(A[3]*D^2)),x,evala); if has(denom(a1),x) or has(denom(a0),x) or degree(a0,x)>2 or degree(a1,x)>2 or (d=3 and coeff(a1,x,2)<>2) or (d=2 and degree(a1,x)>1) then if has([args],`input normalized`) then # This ode is, after normalization, still not an ode with 3 # regular singularities. [] else if d=3 then s:=evala(S+x^2/D); # coefficients of the equation after D -> D+s : v:=[1,2*x^2/D,evala(A[3]-2*S*s+s^2+diff(s,x))] else s:=S; v:=[1,0,C] fi; # multiply by exp(int(s)) to compensate for the D -> D+s : [seq(i*e_int(s,x),i=procname(v,0,C,V,x,`input normalized`))] fi else v:=[seq(seq(coeff(s,x,i),i=0..2),s=[collect(D,x,evala),a0,a1])]; if d=3 then sols3(v,V,x,args) else sols2(v,V,x,args) fi fi end: # 3 finite singularities. Coefficients of equation given in v. sols3:=proc(v,V,x,A) local P,Z,R,S,T; # indicial equation at R. Degree = 2 iff convert(V,`*`) is squarefree. option `Copyright (c) 1997 Waterloo Maple Inc. All rights reserved.`; P:=evala( Z^2*(4*v[2]^3-v[2]^2*v[3]^2-18*v[1]*v[2]*v[3]+4*v[1]*v[3]^3+27*v[1]^2)+(v[2]^2 *v[3]^2-9*v[1]^2-2*v[2]^2*v[3]*R-v[2]*v[3]^2*v[7]-3*v[1]*v[3]*v[7]+R^2*v[2]*v[ 3]*v[8]+16*v[1]*v[2]*v[3]+v[2]*v[3]^2*R*v[8]-4*v[2]^3-4*v[1]*v[3]^3-6*v[2]*v[8 ]*v[1]+2*v[3]^2*v[8]*v[1]-3*v[1]*v[3]*R*v[8]-9*R*v[1]*v[7]+6*R^2*v[2]*v[7]+6*R *v[1]*v[2]+4*v[1]*v[3]^2*R-9*R^2*v[1]*v[8]+12*R^2*v[1]*v[3]-2*R^2*v[3]^2*v[7]-\ 2*v[3]^3*R*v[7]-2*v[2]^2*R*v[8]+7*R*v[2]*v[3]*v[7]-4*v[2]^2*R^2+4*v[2]^2*v[7]) *Z+3*R^2*v[4]-v[3]^2*v[5]*R+2*v[3]*R*v[4]+v[2]*v[5]*R-v[3]*R^2*v[5]+v[2]*v[6]* R^2-3*v[6]*R*v[1]+v[6]*R*v[2]*v[3]-3*v[5]*v[1]+v[6]*v[1]*v[3]-v[3]^2*v[4]+4*v[ 2]*v[4] ); S:=S_and_E(R,Z,P,args,`no indexed RootOfs`); T:=S[2]; # indicates integer difference of exponents. S:=S[1]; S:= [hypergeom([-S[8]+1-S[5]-S[7], S[6]+S[8]+S[4]],[1+S[4]-S[5]],(-x+S[1])*(S[2]-S [3])/(x-S[3])/(S[1]-S[2]))*((x-S[1])/(x-S[3]))^S[4]*((x-S[2])/(x-S[3]))^S[6], ((-x+S[1])*(S[2]-S[3])/(x-S[3])/(S[1]-S[2]))^(-S[4]+S[5])*'hypergeom'([-S[8]+1-S [7]-S[4], S[6]+S[8]+S[5]],[1-S[4]+S[5]],(-x+S[1])*(S[2]-S[3])/(x-S[3])/(S[1]-S [2]))*((x-S[1])/(x-S[3]))^S[4]*((x-S[2])/(x-S[3]))^S[6]]; if T then [S[1],S[1]*_Env_odsolve_int(e_int(-A[2],x)/S[1]^2,x)] else S fi end: sols2:=proc(v,V,x,A) local P,Z,R,S,T,j; # indicial equation at R. Degree = 2 iff convert(V,`*`) is squarefree. option `Copyright (c) 1997 Waterloo Maple Inc. All rights reserved.`; P:=evala( (4*v[1]-v[2]^2)*Z^2+(-2*R*v[7]+2*v[8]*v[1]+v[2]*R*v[8]-4*v[1]-v[7]*v[2]+v[2]^2 )*Z-R*v[5]-v[4]+v[6]*v[1]+v[6]*v[2]*R ); S:=S_and_E(R,Z,P,args,`no indexed RootOfs`); T:=S[2]; # indicates integer difference of exponents. # indicial equation at infinity: P:=evala(Factors(numer(evala(Normal(v[6]+(-v[8]+1)*Z+Z^2))) ,indets([args],RootOf)))[2]; P:=[seq(`if`(has(j,Z),j,NULL),j=P)]; S:=[op(S[1]),evala(RootOf(P[1][1],Z))]; # T:=T or P[1][2]=2 or (nops(P)=2 and type( # evala(Normal(RootOf(P[1][1],Z)-RootOf(P[2][1],Z))),integer)); S:= [hypergeom([-S[5]-S[7]-S[8]+1, S[8]+S[4]+S[6]],[1+S[4]-S[5]],(-x+S[1])/(S[1]-S [2]))*(x-S[1])^S[4]*(x-S[2])^S[6], ((-x+S[1])/(S[1]-S[2]))^(-S[4]+S[5])* 'hypergeom'([S[8]+S[6]+S[5], -S[7]-S[8]+1-S[4]],[1-S[4]+S[5]],(-x+S[1])/(S[1]-S[ 2]))*(x-S[1])^S[4]*(x-S[2])^S[6]]; if T then [S[1],S[1]*_Env_odsolve_int(e_int(-A[2],x)/S[1]^2,x)] else S fi end: S_and_E:=proc(R,Z,P,v,V,x) local j,E,Q,Tr,W,i,intdif,nonintdif; # Now start building the list of singularities: option `Copyright (c) 1997 Waterloo Maple Inc. All rights reserved.`; W:=NULL; for i in V do if degree(i,x)=1 then W:=W,x-i elif degree(i,x)=2 then W:=W,RootOf(i,x),evala(-coeff(i,x,1)-RootOf(i,x)) else if has([args],`no indexed RootOfs`) then # Temporary fix for the situation that the indexed # RootOf's are not yet ready. W:=`DEtools/diffop/split_them`(RootOf(i,x),i) else W:=RootOf(i,x,'index'=1),RootOf(i,x,'index'=2),RootOf(i,x,'index'=3) fi fi od; W:=[W]; E:=NULL; nonintdif:=false; # Now start building the list of exponents: for i in W do Q:=[seq(`if`(has(j,Z),j,NULL),j=evala(Factors( numer(evala(Normal(subs(R=i,P)))),indets([args],RootOf)))[2])]; if Q[1][2]=2 then intdif[i]:=true; Q:=evala(RootOf(Q[1][1],Z)); Q:=Q,Q elif nops(Q)=2 then Q:=evala(RootOf(Q[1][1],Z)),evala(RootOf(Q[2][1],Z)); j:=evala(Normal(Q[1]-Q[2])); if type(j,integer) then intdif[i]:=true; if j<0 then # to make sure that the third entry of # hypergeom gets > 0 Q:=Q[2],Q[1] fi else intdif[i]:=false; nonintdif:=true fi; else intdif[i]:=false; nonintdif:=true; Tr:=-evala(coeff(Q[1][1],Z,1)/coeff(Q[1][1],Z,2)); Q:=RootOf(Q[1][1],Z),Tr-RootOf(Q[1][1],Z) fi; E:=E,Q od; while nonintdif and intdif[W[1]] do W:=[op(W[2..nops(W)]),W[1]]; E:=E[3..nops([E])],E[1],E[2] od; [[op(W),`if`(nops(W)=2,0,NULL),E],not nonintdif] end: # Output: either the singularities, or FAIL if it is not a # Riemann-Papperitz type equation. sing:=proc(C,x) local d,i,v; option `Copyright (c) 1997 Waterloo Maple Inc. All rights reserved.`; d:=denom(C); if degree(d,x)>6 or degree(d,x)-degree(numer(C),x)<2 then # >3 singularities or not regular singular FAIL else v:=[seq(`if`(has(i,x),i,NULL),i=evala(Factors(d,indets([args],RootOf)))[2])]; if max(seq(i[2],i=v))>2 or add(degree(i[1],x),i=v)>3 then # >3 singularities or not regular singular FAIL else # make them monic and return the factors v:=[seq(collect(i[1]/lcoeff(i[1],x),x,evala),i=v)]; while degree(convert(v,`*`),x)<2 do d:=1+max(-1,seq(`if`(type(i-x,integer),i-x,NULL),i=v)); v:=[x-d,op(v)] od; v fi fi end: #savelib('sols','sols2','sols3','S_and_E','sing',\ '`DEtools/RiemannPsols`'):