# Solver for differential equations of order 2 with rational function coeffs # # Finds all solutions that can be expressed in terms of Bessel{I,J,K,Y}(nu, f) # or Whittaker{M,W}(mu, nu, f) or Kummer{M,U}(mu, nu, f) with mu, nu # constants and f an arbitrary rational function. # # Version 2.0.3 --- September 04, 2007 # # Copyright: Mark van Hoeij and Ruben Debeerst. macro( diffop2de = DEtools['diffop2de'], de2diffop = DEtools['de2diffop'], EXP_INT = `DEtools/kovacicsols/e_int`, dchange = PDEtools['dchange'], symmetric_product = DEtools['symmetric_product'], EXP_SOLS = `DEtools/Expsols`, NO_SOL = "No Bessel or Whittaker solutions", KUMMER = `dsolveBessel/KummerScore` ): dsolveBessel := proc() local ans; # use frontend to replace non-algebraic constants by names ans := traperror(frontend(dsolveBessel_doit,[args], [{`+`,`*`,RootOf,radical,nonreal}, indets([args[2..-1]],function)])); if ans=lasterror then if ans=NO_SOL then ans := 0 else error ans fi fi; if _EnvListOutput = true and ans=0 then ans := [] fi; ans end: # Input: # i) a differential operator L and (optional) the domain y # ii) a differential equation L and the dependent variable y # Compute: M, B where M is a first order differential operator and B is a # function of the form C1*BesselI(v,f(x))+C1*BesselK(v,f(x)) such that # M(B) is a solution of L, if such a solution exists. If such a solution M(B) # exists, return it, otherwise return either 0 or NO_SOL. dsolveBessel_doit := proc(L,y) local Dx,x,Sirr,Srest,vflist,vf,ext,T,sqrts,sc,i,j,B,mayB,caseW,caseB; i := NULL; if nargs>1 then if type(y,list(name)) and nops(y)=2 then _Envdiffopdomain:=y elif type(y,function) then x := op(y); _Envdiffopdomain:=[Dx,x]; return procname(de2diffop(L,y)) elif type(y,table) then # This is used to pass along data from kovacicsols T := y[[]]; i := y else error "wrong number or type of arguments" fi fi; if not assigned(_Envdiffopdomain) then error "domain is not assigned" fi; Dx,x := op(_Envdiffopdomain); if degree(L, Dx)<>2 then error "input should be of order 2" fi; ext := indets({args} minus {i}, {RootOf,radical,nonreal}); if ext = {} then Normalizer := normal; else Normalizer := evala; fi; Sirr, Srest, mayB, caseB, caseW := singgenexp(L,ext,T, i); userinfo(1,dsolveBessel,"Generalized exponents ", Sirr,Srest); if mayB then sqrts := SqrtConst(Sirr, ext, T, x); userinfo(1,dsolveBessel,"set of adjusted Sirr and c's",sqrts); vflist := NULL; # vflist contains tuples [v,c,f] for sc in sqrts do i := [findBesselvf(caseB,sc[1],Srest,ext,T)]; vflist := vflist, op(map(x->[x[1],sc[2],x[2]],i)); od; userinfo(1,dsolveBessel,"possible [v,c,f]",vflist); for vf in [vflist] do userinfo(1,dsolveBessel,"testing",vf); B:=besselequiv(L,op(vf)); if B<>0 then return SimplifyAnswer(B) fi od; fi; if caseW>0 then findWhittaker(L, Sirr, Srest, ext, T, x, caseB, caseW) else 0 fi end: # Apply operator exp(Int(d,x)) * B1 on the functions in B2. SimplifyAnswer := proc(d, B1, B2) local c, c2, C, Dx, x, i, de, Y, f, sol, p; global _C1, _C2; userinfo(1,dsolveBessel,"Operator and functions",exp(Int(d,x))*B1, B2); Dx, x := op(_Envdiffopdomain); de := diffop2de(B1,Y(x)); f := eval(de, Y(x) = B2[1]); i := select(has, indets(f,function), op(0,B2[1])); # f := collect(primpart(f, i, 'c'), i, factor); # # Bug, when f := A( 3^(1/2)*(2*x+4) ); i := {f}; # then primpart(f, i, 'c'); returns 1 instead of f. # # Fix: p := proc(f,v) local w,i,j; w := {seq([v[i]=w[i],w[i]=v[i],w[i]], i=1..nops(v))}; w := seq(map(i -> i[j], w), j=1..3); i := primpart(eval(f,w[1]), w[3], 'j'); eval(collect(i,w[3],factor),w[2]), j end: f, c := p(f, i); C := EXP_INT(Normalizer(d + diff(c,x)/c),x); sol := C * f; f := eval(de, Y(x) = B2[2]); i := select(has, indets(f,function), op(0,B2[2])); # f := collect(primpart(f, i, 'c2'), i, factor); # Same bug as above, fix: f, c2 := p(f, i); sol := [sol, Normalizer(c2/c) * C * f]; if _EnvListOutput = true then sol else _C1 * sol[1] + _C2 * sol[2] fi end: # find Whittaker(v,mu,2*f, caseW) findWhittaker := proc(L,Sir,Srest,ext,T,x, caseB, caseW) local vf,v,c,f,Sirr,mult,i,k,s,mulist,ans,mu, vflist, j; if caseW = 3 then for i in Sir do if type(i[3],list) then v := coeff(i[3,1],T,0)/ldegree(i[3,1],T)/2; v := evala(v^2*i[3,2]) else v := coeff(i[3],T,0)/ldegree(i[3],T)/2; v := x^2 - poly_d(v, x, ext) fi; if assigned(c) and evala(c-v)<>0 then error NO_SOL fi; c := v od; Sirr := SqrtConst(Sir, ext, T, x, c)[1,1] else c := 1; Sirr := Sir; fi; vflist := NULL; for i in [findBesselvf(caseB, Sirr,Srest,ext,T)] do if nops(i)=3 then # The logarithmic case for j in i[1..2], [abs(i[1]-1/2), i[2]] do if j[1]=i[3] then vflist := j, vflist else vflist := vflist, j fi od else vflist := vflist, i fi; od; for vf in [vflist] do v, f := op(vf); userinfo(2,dsolveBessel,"Trying v,c,f = ", v,c,f); ans := [1,-1]; for i from 1 to nops(Sirr) while c=1 do mult := -ldegree(Sirr[i,3],T); userinfo(3,dsolveBessel,"[p,mult]",[Sirr[i,1],mult]); k := coeff(Sirr[i,3],T,0); k := [seq(k+j,j=0..mult-1)]; mulist := k/(2*mult); if i=1 then ans := mulist else ans := [seq(`if`(member(true, [seq(seq(type(2*evala(j-k*s),integer) ,k=mulist),s={1,-1})] ),j,NULL),j=ans)] fi; userinfo(4,dsolveBessel,"mu's so far:",ans); if ans = [] then break fi; od; if ans<>[] then f := convert(f,sqrfree,x) fi; userinfo(1,dsolveBessel,"remaining mu's:",ans); for mu in ans do userinfo(1,dsolveBessel,"testing mu",mu); s := sign(length([-mu,-2*f]) - length([mu,2*f])); if s=0 then s := `if`(abs_value(mu)=mu, 1, -1) fi; i := whittakerequiv(L, s*mu, v, s*2*f, c); if i<>0 then return SimplifyAnswer(i) fi od; od; 0 end: whittakerequiv := proc(L,mu,v,f,c) local Dx,x,eq,t,L2,Y,g,wc, s,w; global KUMMER; Dx,x := op(_Envdiffopdomain); if c=1 then # See if the answer is shorter in terms of Kummer functions: g := seq(seq([1/2-s*mu+w, 1+2*w,s*f] ,w=[v,-v]), s=[1,-1]); g := sort([g],length)[1]; if length(g) < length([mu,v,f]) + 5 * KUMMER then return Kummerequiv(L,op(g),x) fi fi; # Whittaker equation after subs: (x,mu) -> sqrt(c) * (x,mu) eq :=diff(Y(x),x,x)-(4*evala(v^2)+c*x^2-4*c*mu*x-1)/4/x^2*Y(x); eq := dchange({x=subs(x=t,f)},eq,[t], params=(indets([args],{name,symbol}) minus {x}) ); L2 := de2diffop(eq,Y(t)): L2 := equiv(L2,L); if L2 <> 0 then wc := sqrt(c, 'symbolic'); t:=1/sqrt(1/c,'symbolic');if length(t) 0 then L2, [KummerM(mu, nu, f), KummerU(mu, nu, f)] else 0 fi end: # Input: differential operator L, extensions ext, name T # Output: a list of singularities S such that # for each singularity s of L at p, S contains [p,t,d,pol,deg], where # - t is local parameter, ie x-p or 1/x # - if D is exponent difference in terms of T, then d is # d=D if D\in K(p,t) # d=[c,d'] if D=c*sqrt(d') and has(c,T) # d=x^2-D if exp difference is sqrt(D)\nin K(p) # d=pol if D\nin K, pol is polynomial with zero D # - pol is polynomial over K with zero p # - deg = degree(pol) singgenexp := proc(L,ext,T, AA) local Dx,x,p,S,Sirr,Srest,g,t,d,c,q,R,extp,mayB,caseB,caseW,muD,Sr,A; global KUMMER; # Used to see if we should use Kummer functions. mayB := true; # true means L may have Bessel solutions Sirr := NULL; Srest := NULL; Sr := NULL; caseB, caseW, muD, KUMMER := {}, {}, 0, 0; Dx,x := op(_Envdiffopdomain); A := evala(Primpart(L,Dx)); A := [seq(coeff(A,Dx,i), i=0..2)]; p := A[3]; userinfo(3,dsolveBessel,"singularities are zeros of",p); if nargs=3 then # Each finite non-apparant singularity is a root of this p: # # p := evala(Gcd(p, denom(Normalizer( # DEtools[LCLM](L,Dx-rand(10..100)()) )) )); # # The following line computes such p much faster: # p := seq(`if`(A[2]=0 or q[2]>1, q[1], Normalizer(q[1] / gcd( gcd(q[1], A[2] + diff(A[3],x)), diff(A[1],x)*A[2]-A[1]*diff(A[2],x)-A[1]^2) )), q=sqrfree(A[3],x)[2]); p := mul(q, q=[p]); S := [infinity, seq(i[1],i=evala(Factors(p,ext))[2])] else # S was passed along, no need to recompute it: S := AA[{}] fi; # Every non-apparant singularity is in S for q in S do if q=infinity then p := infinity else if not has(q,x) then next fi; p := RootOf(q,x) fi; userinfo(3,dsolveBessel,"processing singularity",p); extp := indets(p,RootOf) union ext; if nargs=4 and assigned(AA[p]) then # This info was passed along, no need to recompute it: g := [AA[p]] else # g := DEtools[gen_exp](L,[Dx,x],T,x=p); # # The above line works OK, however, kovacicsols has # a gen_exp implementation optimized for order 2. # So we'll use that instead: # g := [`DEtools/kovacicsols/gen_exp`(A,T,x,p,extp)] fi; userinfo(3,dsolveBessel,"point and gen_exp",[p,g]); t := `if`(p=infinity, 1/x, x-p); # t = local parameter at p if nops(g)=1 and nops(g[1]) = 2 then # gen_exp involves an algebraic extension if lhs(g[1,-1])<>T then error NO_SOL fi; R := indets(g,RootOf) minus extp; if nops(R)<>1 then error "should not happen" fi; R := R[1]; c := coeff(collect(g[1,1],R),R); d := evala(discrim( evala(Norm(x-R,extp,extp)), x)); if has(c,T) then d := [collect(c,T,evala),d] else d := x^2-evala(c^2*d); fi else d := g[1,1] - `if`(nops(g)=2, g[2,1], g[1,2]); d := collect(d, T, evala); if not has(d,T) and indets(d,RootOf) minus ext<>{} then d := poly_d(d, x, ext) else d := abs_value(d) fi; fi; if d=0 or ( type(d,integer) and d<>0 and DEtools['formal_sol'](L,`has logarithm?`,x=p) ) then # d is an integer and has a logarithm iff # there is a gauge operation such that d=0 d := {d} fi; p := [p,t,d,op(`if`(q=infinity, [1,1], [q,degree(q,x)]))]; userinfo(2,dsolveBessel,"gen_exp-array", p); if has(d,T) then if type(d,list) then # d in algebraic extension t := coeff(d[1],T,0); if t=0 then # No Whittaker solutions: caseW := {0,1} else mayB := false; # 0: no Whittaker solutions # 1: when mu is rational # 2: when mu is irrational but in K # 3: when mu is not in K caseW := caseW union {3} fi else t := coeff(d,T,0); mayB := mayB and type(t,integer); if type(t,rational) then muD := igcd(muD, denom(t)*ldegree(d,T)); # denom(2*mu) divides muD, so if muD=1 # then 2*mu\in Z, which is Bessel case caseW := caseW union {1,`if`(muD=1,0,1)} elif indets(t,RootOf) minus ext = {} then caseW := caseW union {2} else caseW := caseW union {3} fi; # if generalized exponents look like Kummer # after change of variable, then increase the # Kummer score; if it looks like whittaker # then decrease the Kummer score: c := {seq(min(0,ldegree(t[1],T)),t=g)}; KUMMER := KUMMER + p[5]* `if`(nops(c)=1, op(c), abs(c[1]-c[2])) fi; if not mayB then userinfo(1,dsolveBessel, "There are no Bessel solutions") fi; if nops(caseW)>1 then userinfo(1,dsolveBessel, "There are no Whittaker solutions"); fi; Sirr := Sirr, p else if type(d,set) then #log case caseB := caseB union {0} elif has(d,x) then caseB := caseB union {3} elif not type(d,rational) then caseB := caseB union {2} elif not type(d,integer) then caseB := caseB union {1} else if d>2 then # Apparant singularities with large d Sr := Sr,p fi; next fi; Srest := Srest, p; fi; if nops(caseB)>1 or (nops(caseW)>1 and not mayB) then error NO_SOL fi; od; if Sirr = NULL then error NO_SOL elif caseB={} then Srest := Sr fi; [Sirr], [Srest], mayB, caseB, `if`(nops(caseW)>1, 0, op(caseW)) end: # Minpoly of c should be an integer shift of x^2-const, if so, return x^2-const poly_d := proc(c, x, ext) local t,d; d := sqrfree(evala(Norm(x-c,ext,ext)),x)[2,1,1]; d := collect(d/lcoeff(d,x),x,Normalizer); t := coeff(d,x,1)/2; if degree(d,x)<>2 or not type(t,integer) then error NO_SOL fi; collect(subs(x=x-t,d),x,Normalizer) end: # Return -a if a looks negative, otherwise return a. abs_value := proc(a) local i; i := traperror(sign(a)); if i=lasterror then i := traperror(sign(Re(evalf[10](a)))) fi; if i=-1 then -a else a fi end: # Input: The list of irregular singularities Sirr. # Let K denote the field given by ext, and let T be the name used in # gen_exp, and x the independent variable. # # Each i in Sirr is a list with # i[1] = point p in C union {infinity} # i[4] = minpoly of i[1] over K, whenever i[1] is not infinity. # i[5] = degree of this minpoly # i[2] = local parameter (not used in this procedure) # i[3] = Either the gen.exp-difference D in K(p)[1/T], # or a list [E,D] such that gen.exp-difference is E * sqrt(D) # and where E in K(p)[1/T] and D in K(p). # # Output: a set of pairs [Updated(Sirr), c]. # Here an updated Sirr is one in which i[3] is never a list [E,D], but # always an element D in K(p)[1/T], and where the gen.exp-difference i[3] # (which was E*sqrt(D) or D depending on whether i[3] was a list or not) # has been divided by sqrt(c) for all i in Sirr. # This means that if we had an i[3] = [E,D] in the input, then for # this D it must have been true that D/c is a square in K(p). And if we # had an i[3] = D in the input, then for this i, we must have had sqrt(c) # in K(p) where p denotes i[1]. # # In most cases the c can be read off from one of the i[3]'s. In the # remaining case we'll have to find candidate c's by computing Subfields. # # The c's coming out of this procedure are used as follows: When we compute # the f for which the solution is in terms of BesselI(v, f) then this f # should still be multiplied by sqrt(c). If this c is negative # then replace c by -c and replace BesselI by BesselJ. # SqrtConst := proc(Sirr, ext, T, x) local Sodd, i, nl, ol, S, disc, ans, c, flag, res, extp, pol, wc, i3; # The [Updated(Sirr), c]'s found so far: ans := NULL; if nargs=5 then # c is already given in the optional input args[5] S := {args[5]} else # Singularities for which K(p) has odd degree over K: Sodd := [seq(`if`(irem(i[5],2)=1, i, NULL), i=Sirr)]; # The number of i's for which i[3] is a list: nl := add(`if`(type(i[3],list),1,0),i=Sirr); ol := add(`if`(type(i[3],list),1,0),i=Sodd); # c can only be 1 if no i[3] is a list (i.e. if nl=0). if nl=0 then ans := [Sirr, 1] fi; # If any odd sing does not have a list, then no cases other than # perhaps c = 1 need to be considered: if ol < nops(Sodd) then return {ans} fi; # At this point i[3] must be a list for all i in Sodd. # Find all singularities where c is easy to read from: S := NULL; for i in Sirr do if not type(i[3],list) then # Such i is not in Sodd, so the degree is even. If # the degree is 2, then c is easy to read. if i[5]=2 then S := S, Normalizer(discrim(i[4],x)) fi; next fi; disc := i[3][2]; if indets(disc,RootOf) minus ext <> {} then if member(i, Sodd) then disc := evala(Norm(disc, ext, ext)) else next fi fi; S := S, Normalizer(disc) od; #if S = NULL then # We found no points where c was easy to read, this # means we'll have to resort to Subfields to get the # potential values of c. #S := seq(`if`(type(i[3],list), # sqrfree(evala(Norm(x^2-i[3][2],ext,ext)),x)[2,1,1], i[4]) #,i=Sirr); S := seq(`if`(type(i[3],list), subs(_Z=x,op(1, lhs( evala(Primfield({`if`(degree(i[4],x)>0, RootOf(i[4],x), NULL), RootOf(sqrfree(evala(Norm(x^2-i[3][2],ext,ext)),x)[2,1,1],x) }))[1,1] ))) , i[4]) ,i=Sirr); S := evala(Subfields({S}, 2, ext, x)); S := {seq(evala(Norm(x-i,ext,ext)), i=S)}; userinfo(2,dsolveBessel,"subfields",S); #else # # Found at least some point where c could be read, so # # just pick the first one and try to simplify it a bit: # S := sort([S],length)[1]; # S := {PolynomialTools:-Shorten(x^2 - S,x)} #fi; # The set of all possible c's other than c=1. S := {seq(Normalizer(discrim(i,x)/4), i=S)}; userinfo(2,dsolveBessel,"possible c's other than c=1 in f=sqrt(c)*g",S); fi; for c in S do # adjust Sirr if possible flag := true; res := NULL; for i in Sirr do extp := ext union indets(i[1],RootOf); pol := `if`(type(i[3],list), i[3,2], 1)*x^2 - c; pol := evala(Primpart(pol,x)); pol := evala(Factors(pol, extp))[2]; if nops(pol)=1 then flag := false; break fi; pol := collect(pol[1,1],x); wc := coeff(pol,x,0)/coeff(pol,x,1); i3 := `if`(type(i[3],list), i[3,1], i[3]) / wc; res := res, subsop(3 = collect(i3,T,evala), i) od; if flag then ans := ans, [[res], c] fi od; {ans} end: # Input: list of singularities Sirr and Srest, extensions ext, name T # Output: list of possible pairs [v,f] findBesselvf := proc(caseB, Sirr, Srest, ext, T) local f, inf; f := besselsubst(Sirr,T,ext); userinfo(1,dsolveBessel,"Possibilities for f ",f); inf := member(infinity, map(i->i[1],Sirr)); if has({2,3},caseB) then # find bessel function with v \nin \Q findBesselvfirrat(Srest,f,inf,ext) elif caseB={0} then # find bessel function with v=0 findBesselvfln(Srest,f,ext) elif caseB={1} then # find bessel function with v\in\Q findBesselvfrat(Srest,f,ext) else # find bessel function with v\in\Q where all singularity # differences are integers findBesselvfint(f,inf,ext,Srest) fi; end: # Input: set of singularities SC with exp differences not in K, # set of possible rational functions f, a boolean inf=(\inf\in\Sirr), # extensions ext # Output: list of possible pairs [v,f] with v\nin K # such that exp-gauge-transformations of B(v,f(x)) can be a solution # of the differential operator belonging to SC findBesselvfirrat := proc(SC,f,inf,ext) local Dx,x,z,i,g,gg,h,hh,d,r,a,b,l,v,sol,vf,ordeRoot,hasx; Dx,x := op(_Envdiffopdomain); vf := NULL; #gen-exp differences are sqrt(.)+rational #make minpolys for sqrt(.) hasx := 0; for i from 1 to nops(SC) do d[i] := SC[i,3]; if has(d[i],x) then hasx := hasx+1 fi; od; if hasx = 0 then return findBesselvfK(args) elif hasx < nops(SC) then error NO_SOL fi; userinfo(3,dsolveBessel,"minpolys for gen-exp differences", [seq(d[i],i=1..nops(SC))] ); # all differences are a rational multiple of d[1] # compute these rational coefficients a := coeff(d[1],x,0); r[1] := 1; for i from 2 to nops(SC) do b := coeff(d[i],x,0); r[i] := sqrt(Normalizer(b/a)); if not type(r[i],rational) then return NULL fi; od; userinfo(4,dsolveBessel,"ratios for gen-exp differences", [seq(r[i],i=1..nops(SC))] ); l := ilcm( seq(denom(r[i]),i=1..nops(SC)) ); h := mul(SC[i,4]^(r[i]*l),i=1..nops(SC)); userinfo(2,dsolveBessel,"numerator must be a power of",h); # test all possibilities of f (try +/- each polar part) for i from 1 to 2^(nops(f)-1) do g := Normalizer(possibility(f,i)); userinfo(3,dsolveBessel,"Testing possibility", g); # get c from first singularity in SC g := changeconstant(g,x,SC[1,1],ext); if g=NULL then next; fi; # test if other singularities in SC are zeros if not testzeros(g,x,map(j->j[1],SC)) then next; fi; gg:=numer(g); if degree(h,x)=0 then a:=1; else a:=degree(gg,x)/degree(h,x); fi; # now h^a=g if not type(a,integer) then next fi; gg:=collect(gg/lcoeff(gg,x),x,Normalizer); hh:=collect(h/lcoeff(h,x),x,Normalizer); if Normalizer(gg-hh^a)=0 then if SC[1,1]=infinity then ordeRoot := degree(denom(g),x)-degree(gg,x) else # h^a contains i'th factor a*r[i]*l times, and r[1]=1 ordeRoot := a*l fi; v := sqrt(factor(-coeff(d[1],x,0)), 'symbolic'); v := v / (2*ordeRoot); vf := vf, [v,g]; fi; od; vf end: # Input: set of singularities SC with non-rational exp differences in K, # set of possible rational functions f, a boolean inf=(\inf\in\Sirr), # extensions ext # Output: list of possible pairs [v,f] # such that exp-gauge-transformations of B(v,f(x)) can be a solution # of the differential operator belonging to SC findBesselvfK := proc(SC, f, inf, ext) local Dx,x,z,i,g,gg,h,hh,r0,r1,a,l,v,vv,k,sol,vf,ordeRoot,flag; Dx,x := op(_Envdiffopdomain); vf := NULL; # all differences are a rational multiple of d[1] # compute these rational coefficients r0[1] := 0; r1[1] := 1; for i from 2 to nops(SC) do r0[i],r1[i] := compare(SC[i,3],SC[1,3]); od; userinfo(4,dsolveBessel,"ratios for exp differences", [seq(r1[i],i=1..nops(SC))] ); l := ilcm( seq(denom(r1[i]),i=1..nops(SC)) ); h := mul(SC[i,4]^(r1[i]*l),i=1..nops(SC)); userinfo(2,dsolveBessel,"numerator must be a power of",h); # test all possibilities of f for i from 1 to 2^(nops(f)-1) do g := Normalizer(possibility(f,i)); userinfo(3,dsolveBessel,"Testing possibility", g); # get c from first singularity in SC g := changeconstant(g,x,SC[1,1],ext); if g=NULL then next; fi; # test if other singularities in SC are zeros if not testzeros(g,x,map(j->j[1],SC)) then next; fi; gg:=numer(g); if degree(h,x)=0 then a:=1; else a:=degree(gg,x)/degree(h,x); fi; if not type(a,integer) then next fi; gg:=collect(gg/lcoeff(gg,x),x,Normalizer); hh:=collect(h/lcoeff(h,x),x,Normalizer); if Normalizer(gg-hh^a) <> 0 then next fi; # now h^a=gg if SC[1,1]=infinity then ordeRoot := degree(denom(g),x)-degree(gg,x) else # h^a contains the i'th factor a*l times ordeRoot := a*l fi; v:=SC[1,3]/(2*ordeRoot); vv := seq(v+k/(2*ordeRoot),k=0..(2*ordeRoot-1)); userinfo(2,dsolveBessel,"set for v", vv); # test v with other zeros for v in vv do flag := true; for k from 2 to nops(SC) do ordeRoot := `if`(SC[k,1]=infinity, degree(denom(g),x) - degree(gg,x), r1[k]*a*l); if not type(Normalizer(v*2*ordeRoot-SC[k,3]) ,integer) then flag := false; break fi; od; if flag then vf := vf, [v,g] fi; od; od; vf end: # Output: a new variable newvar := proc() local t; t end: # Input: A,B in k # Output: r0,r1 such that A=r0+r1*B compare := proc(A,B) local iszero,E,sol; iszero := A- (r0+r1*B); iszero := evala(Expand(numer(Normalizer(iszero)))); iszero := convert(iszero,RootOf); do E := indets(iszero,RootOf); if E = {} then break fi; iszero := subs(E[1]=newvar(),iszero); od; sol := SolveTools:-Linear({coeffs(expand(iszero), indets(iszero,{name,symbol}) minus {r0,r1})},{r0,r1}); if not type(sol,set) or not type(map(rhs,sol),set(rational)) then error NO_SOL fi; op(subs(sol, [r0,r1])) end: # Input: set of singularities Sln with exp differences in Z and log in # formal solution, set of possible rational functions f, a boolean # inf=(\inf\in\Sirr), extensions ext # Output: list of possible pairs [v,f] with v=0 # such that exp-gauge-transformations of B(v,f(x)) can be a solution # of the differential operator belonging to SC findBesselvfln := proc(Sln,f,ext) local Dx,x,i,j,g,sol,v,vf; Dx,x := op(_Envdiffopdomain); vf := NULL; # test all possibilities of f for i from 1 to 2^(nops(f)-1) do g := Normalizer(possibility(f,i)); userinfo(3,dsolveBessel,"Testing possibility", g); # compute constant g := changeconstant(g,x,Sln[1,1],ext); if g=NULL then next; fi; # test for other zeros if not testzeros(g,x,map(j->j[1],Sln)) then next; fi; v := 2*max(degree(numer(g),x),degree(denom(g),x)); v := add(j[3,1]*j[5], j=Sln)/v; vf := vf, [ceil(v),g,v] od; vf end: # Input: set of singularities SQ with exp differences in Q, # set of possible rational functions f, extensions ext # Output: list of possible pairs [v,f] with v in Q # such that exp-gauge-transformations of B(v,f(x)) can be a solution # of the differential operator belonging to SC findBesselvfrat := proc(SQ,f,ext) local Dx,x,i,j,k,g,h,fctrs,mult,v,vv,c,sol, newh,vf; Dx,x := op(_Envdiffopdomain); vf := NULL; for i from 1 to 2^(nops(f)-1) do g := possibility(f,i); userinfo(3,dsolveBessel,"Testing possibility", g); # get constant from first singularity g := changeconstant(g,x,SQ[1,1],ext); if g=NULL then next; fi; if not testzeros(g,x,map(j->j[1],SQ)) then next; fi; h := numer(g); for j from 1 to nops(SQ) do # look for multiplicities of the zeros mult[j] := 0; if SQ[j,1] = infinity then mult[j] := degree(denom(g))-degree(numer(g)); else mult[j] := 0; while evala(Rem(h, SQ[j,4], x, 'newh')) = 0 do mult[j] := mult[j] + 1; h := newh od; # mult[j] \neq 0 since we already know # that this is a zero fi; # possibilities for v v := [seq(op([SQ[j,3]-i,SQ[j,3]+i]),i=0..mult[j])] / (2*mult[j]); if j=1 then vv := [seq(`if`(MatchInt([0,1/2,op(v[1..i-1])], v[i]), NULL, abs(v[i])), i=1..nops(v))] else vv := [seq(`if`(MatchInt(v,i),i,NULL),i=vv)] fi; if nops(vv) = 0 then break; fi; od; userinfo(3,dsolveBessel,"multiplicities of factors", [seq(mult[j],j=1..nops(SQ))]); userinfo(2,dsolveBessel, "possibilities for v", vv); # if we get here, then we have a function g and a parameter v # where for each singularity [x,t,d] in SQ there is a # corresponding factor t in g, and it's multiplicity m is # such that d=2*v*m. The rest-polynomial h has now to be a # p-power such that 2*v*p \in \Z for v in vv do # rest must be a p-power c:=lcoeff(h,x); newh:=ispower(h/c,x,denom(2*v)); if Normalizer(h/c-newh^denom(2*v)) = 0 then vf := vf, [v,g]; fi; od; od; vf end: MatchInt := proc(L, i) local j; for j in L do if type(i-j,integer) or type(i+j,integer) then return true fi od; false end: # Input: set of possible rational functions f, a boolean inf=(\inf\in\Sirr), # extensions ext # Output: list of possible pairs [v,f] with v in Q # such that exp-gauge-transformations of B(v,f(x)) can be a solution # of a differential operator with no half-regular singularities findBesselvfint := proc(f,inf,ext, Sr) local Dx,x,ff,n,p,i,g,lc,a,k,v,sol,C,c,cc,fctrs,vf,s; Dx,x := op(_Envdiffopdomain); vf:= NULL; for i from 1 to 2^(nops(f)-1) do ff := Normalizer(possibility(f,i)+C); userinfo(3,dsolveBessel,"Testing possibility", ff); if inf then n := degree(numer(ff)); else n := degree(denom(ff)); fi; for p from 2 to n do if irem(n,p)=0 then userinfo(3,dsolveBessel, "Test whether numerator is a p-power", p); # if inf then test p-power and get c # if not then get c from p power or c is zero g := numer(ff); lc := lcoeff(g,x); a := ispower(collect(g/lc,x,Normalizer),x,p); userinfo(5,dsolveBessel, "Is a power of", a); s := numer(Normalizer(evala(Content(g-lc*a^p,x)))); s := evala(Factors(s,ext)); cc := {seq(`if`(degree(k[1],C)=1, -coeff(k[1],C,0)/coeff(k[1],C,1),NULL),k=s[2])}; if not inf then # test c=0 fctrs := sqrfree(subs(C=0,g),x); if add(modp(k[2],p), k=fctrs[2])=0 then # the multiplicity of each factor of g # is a multiple of p cc:=cc union {0}; fi; fi; userinfo(2,dsolveBessel, "possibilities for c", cc); for c in cc do g := subs(C=c,ff); for k from 1 to p-1 do v := k/(2*p); if denom(2*v)

[] and UpdateV(v,g,x,Sr,'v') then vf := [v,g], vf else vf := vf, [v,g] fi; od; od; fi od; od; vf end: # Update v in such a way to make sure that no gauge transformation # is used when none is needed. This procedure is not needed for # correctness or completeness, but only to make the output look nicer. UpdateV := proc(v,g,x,Sr,newv) local s,m,mg,h,N,n1,n2; m := {op(map(s -> s[3], Sr))}; m := min(op(m)); for s in Sr while s[3]>m do od; if s[1]=infinity then mg := degree(denom(g),x)-degree(numer(g),x) else mg := 0; h := numer(g); while evala(Rem(h, s[4], x, 'h'))=0 do mg := mg+1 od fi; if mg=0 then return false fi; n1 := solve( (v+N)*2*mg = m, N); if type(n1,integer) then newv := v+n1; return true fi; n2 := solve( (N-v)*2*mg = m, N); if type(n2,integer) then newv := n2-v; return true fi; if n1>1 then if abs(round(n1)-n1) <= abs(round(n2)-n2) then newv := v+round(n1) else newv := round(n2)-v fi fi; false end: # Input: set of singularities Sirr of a differential operator L with # non-constant generalized exponent # Output: the set {f1,...fn} such that if L has a Bessel solution, then it is # of the form r1*E(x)*B(v,f(x))+r2*(E(x)*B(v,f(x)))', where # f=(-1)^k1*f1+...+(-1)^kn*fn+C, r1,r2\in\C(x), ki\in\Z and E(x) is a # solution of a first order differential equation besselsubst := proc(Sirr,T,ext) local i,d,f,ff,j; f := NULL; for i from 1 to nops(Sirr) do d := 1/2 * Sirr[i,3]; ff := add( coeff(d,T,j)/j*Sirr[i,2]^j, j=ldegree(d,T) .. -1); f := f, evala(Trace(ff,ext,ext)) od; [f] end: # Input: a list f=[f1,...fn] and an integer 1<=k<=2^n # Output: the sum +-f1+-...+-fn, where the k-th possibility # of signs is choosen possibility := proc(f,k) local b,n,a,i; n := nops(f); if k>2^n then error "there are just 2^n possibilities" fi; b := k-1; for i from n by -1 to 1 do b := iquo(b,2,a[i]) od; add((-1)^a[i]*f[i],i=1..n) end: # Input: A monic polynomial f\in K[x] and a positive integer p\in\N # Output: g\in K[x], s.t. # if there exists a solution to y^p=f, then y=g is such a solution ispower := proc(f,x,p) local d,n,A,Ap,a,i,b; if p=1 or degree(f,x)=0 then return f fi; d:=degree(f,x); n:=d/p; if not type(n,integer) then return NULL fi; A:=add(a[i]*x^i,i=0..n-1)+x^n; Ap:=collect(A^p,x); for i from 1 to n do # Solve a linear equation in one unknown: b[n-i]:=solve(subs({seq(a[n-j]=b[n-j],j=1..i-1)},coeff(Ap,x,d-i)) -coeff(f,x,d-i),a[n-i] ); od; Ap:=subs({seq(a[n-j]=b[n-j],j=1..n)},Ap); subs({seq(a[n-j]=b[n-j],j=1..i-1)},A) end: # Input: a function f(x) and a point p # Output: f(x)+c such that f(p)+c=0, # if p=\infty then f(x) is returned changeconstant := proc(f,x,p,ext) local c,pp,ff,C; if p=infinity then return f; fi; if type(p, RootOf) then pp:=evala(Norm(x-p,ext,ext)); ff:=numer(Normalizer(f+C)); c:=SolveTools:-Linear({coeffs(evala(Rem(ff,pp,x)),x)},{C}); if c=NULL then NULL else Normalizer(f+subs(c[1],C)) fi; else Normalizer(f-subs(x=p,f)) fi; end: # Input: a function f(x) and a set of points pp # Output: true, if all points in pp are a zero of the numerator # and false otherwise # for \infty the condition changes to # degree(denom(f))>degree(numer(f)) testzeros := proc(f,x,pp) local p,g; g:=numer(f); for p in pp do if p=infinity then if not degree(g,x) sqrt(c) * x eq := x^2*diff(Y(x),x,x) + x*diff(Y(x),x) - (c*x^2 + evala(v^2))*Y(x); cpos := abs_value(c); modifiedbessel := evalb(c = cpos); eq := dchange({x=subs(x=t,f)},eq,[t], params=(indets([args],{name,symbol}) minus {x})); L2 := de2diffop(eq,Y(t)): L2 := equiv(L2,L); if L2 <> 0 then if type(v,rational) and nargs=4 and L2[2]<>1 and denom(v)>2 then g := length(L2[2]); t := procname(L, abs(v+1-2*frac(v)), c, f, g); if length(t[2])2 or degree(M2,Dx)<>2 then error "case not handled yet" fi; T1:=Normalizer(coeff(M1,Dx,1)/coeff(M1,Dx,2)/2); L1:=symmetric_product(M1,Dx-T1); T2:=Normalizer(coeff(M2,Dx,1)/coeff(M2,Dx,2)/2); L2:=symmetric_product(M2,Dx-T2); C0:=Normalizer(coeff(L1,Dx,0)/coeff(L1,Dx,2)); D0:=Normalizer(coeff(L2,Dx,0)/coeff(L2,Dx,2)); if not Same_p_curvature(C0,D0,x) then userinfo(2,dsolveBessel,"case discarded by p-curvature test"); return 0 fi; D0:=Normalizer(D0-C0); if D0=0 then userinfo(1,dsolveBessel,"No gauge transformation needed"); d := 0; S := 1 else userinfo(3,dsolveBessel,"Computing gauge transformation"); C1:=Normalizer(diff(C0,x)); D1:=Normalizer(diff(D0,x)); d:=Normalizer(D1/D0); # v = coeff_list of symmetric_product(L1, L2) v:=map(Normalizer,[2*diff(C1,x)+diff(D1,x)-2*d*C1-D1*d+D0^2, 6*C1+D1-4*d*C0, 2*(D0+2*C0), -d, 1]); v, d := RemoveAppSingOrd4(v,x); # This _Env passes information on apparant singularities # to the expsols code, which then runs faster: _Env_DF_sq := table([x, denom(C0)*denom(D0)]); # The options tell Expsols to only look at exponential # solutions with generalized exponents in 1/2 * Z (i.e. only # exponential solutions whose square is a rational function) v:=EXP_SOLS(v,0,x,`use Int`,`no algext`,radical,denom,2); if nops(v)=0 then return 0 elif nops(v)>1 then userinfo(1,'dsolve',"input is reducible"); return 0 fi; v1 := `if`(type(v[1],`*`),[op(v[1])],v); v := 1; for i in v1 do if not has(i,x) then next fi; if type(i,function) and op(0,i)=exp and type(op(i),function) and op(0,op(i))=Int and op(2,op(i))=x then d := d + op(1,op(i)) else v := v * i fi od; w[0] := Normalizer(v); for i to 3 do w[i] := Normalizer(diff(w[i-1],x) + d*w[i-1]) od; S:=collect(w[0]*Dx + (w[3]+(4*C0+5*D0)*w[1]+w[0]*(2*C1+D1))/2/D0-2*w[1], Dx,Normalizer); S := symmetric_product(S, Dx+T1); S := evala(Primpart(S, Dx, 'i')); d := d + diff(i,x)/i; S := frontend(collect,[S,Dx,factor]) fi; d := Normalizer(d + T1 - T2); d, S # This represents the operator exp(Int(d,x)) * S end: # Try to remove some apparant singularities we may have introduced. RemoveAppSingOrd4 := proc(v, x) local d,i,Dx,L; d := mul(`if`(i[2]=1,i[1],1),i=sqrfree(denom(v[-2]),x)[2]); if not has(d,x) then return v, 0 fi; i := denom(Normalizer(v[-2]-2*diff(d,x)/d)); d := denom(Normalizer(i/d)); if not has(d,x) then return v, 0 fi; i := 2*v[-3] - 3*(v[-2]*diff(d,x)+2*diff(d,x,x))/d + (3*diff(d,x)/d)^2; d := denom(Normalizer( denom(Normalizer(i))/d )); if not has(d,x) then return v, 0 fi; d := 1/2 * diff(d,x)/d; L := symmetric_product(add(v[i+1]*Dx^i,i=0..4), Dx-d, [Dx,x]); [seq(Normalizer(coeff(L,Dx,i)),i=0..4)], -d end: # Operators Dx^2 + a, Dx^2 + b are compared mod 3 # If the reduction mod 3 fails, then compare mod 5 instead. Same_p_curvature := proc(a,b, x) local v,i,s,s3,T; if a=b then return true fi; v := traperror(`mod/ReduceField`([a,b],{x},3)); if v = lasterror then return Same_5_curvature(args) fi; for i to 2 do s := Normal(v[i]) mod 3; s3 := subs({seq(j=j^3, j=indets(s, {name,symbol,RootOf}))}, s); s := Normal(diff(s^2,x,x)) mod 3; T[i] := s3 + s od; evalb(Normal(T[1] - T[2]) mod 3 = 0) end: # Operators Dx^2 + a, Dx^2 + b are compared mod 5 Same_5_curvature := proc(a,b, x) local v,i,s,s5,j,T; v := traperror(`mod/ReduceField`([a,b],{x},5)); if v = lasterror then return true fi; for i to 2 do s := Normal(v[i]) mod 5; s5 := subs({seq(j=j^5, j=indets(s, {name,symbol,RootOf}))}, s); s := Normal( s * (2*s^2 + diff(s,x,x)) ) mod 5; for j to 4 do s := Normal( diff(s,x) ) mod 5 od; T[i] := s + s5 od; evalb(Normal(T[1] - T[2]) mod 5 = 0) end: