with(DEtools): _Envdiffopdomain := [Dx,x]: P2 := x*(x-1)*(64*x-1)*(3*x-2)*(6*x+1)*Dx^2 + (4608*x^4 - 6372*x^3 + 813*x^2 + 514*x-4)*Dx + 4*(576*x^3 - 801*x^2 -108*x + 74); singularities := {solve(lcoeff(P2,Dx))} union {infinity}; for i in singularities do lprint(i); print( gen_exp(P2,T,x=i) ); print( formal_sol(P2,`has logarithm?`,x=i) ); od: # Shows that all exponent-differences are integers, and that x = -1/6 is a false singularity. TrueSing := singularities minus {-1/6}; # In this file we compute the (0,1/3,0) type solutions. # Then all roots and all poles of f must be in the set TrueSing, because all # singularities of P2 are logarithmic. # Find the products of x-i with i in S, with degree N. FindProducts := proc(S, N, x) FindProducts_Num_Deg(S, N, N, x) end: # Find all products with degree(numerator)=degN and degree(denominator)=degD # such that every element of S appears either in the denominator or in the numerator FindProducts_Num_Deg := proc(S, degN, degD, x) local s,T,t; options remember; if S={} then if degN=0 and degD=0 then {1} else {} fi else s := S[1]; T := S minus {s}; if s=infinity then t:=1 else t:=x-s fi; {seq(seq(t^j * k, k=procname(T, degN-j, degD, x)), j=1..degN), seq(seq(t^j * k, k=procname(T, degN, degD+j, x)), j=-degD..-1)} # replacing that last -1 by a 0 removes the condition that every element of S # must appear in the denominator (that's used for (0,0,0) exponent-differences). fi end: IsCube := proc(f,x,unknown) local d,u,q,S; d := degree(f,x)/3; q := add(u[i]*x^i, i=0..d-1)+x^d; S := solve({coeffs(collect(rem(f,q^3,x),x),x)}, {unknown} union {seq(u[i],i=0..d-1)}); if S={} or S=NULL then NULL else eval(unknown, S) fi end: # Search for (0,1/3,0) type solutions. found := false; for degf in [3,6,9,12,15,18,21,24] # if degf>24 we give up while not found do PR := FindProducts(TrueSing,degf,x); for i in PR do f := c*i; # For (0,1/3,0) we need the numerator of f-1 to be a cube v := IsCube(numer(f-1),x,c); if v<>NULL then lprint(factor(eval(f, c=v))); found := true fi od: od: # We found: f := -27*x*(3*x-2)/(x-1)^2/(64*x-1); # (we also found its reciprocal, that's because (0,1/3,0) has a symmetry x -> 1/x). Le := Dx^2-(-2*x+x*e0+x*e1+1-e0)/x/(x-1)*Dx+1/4*(e0-1+e1+ei)*(e0-1+e1-ei)/x/(x-1); # Has exponent-differences (e0,e1,ei). One of its solutions is: sol1 := hypergeom([1/2-1/2*e0-1/2*e1+1/2*ei, 1/2-1/2*e0-1/2*e1-1/2*ei],[1-e0],x); L030 := subs(e0=0, e1=1/3, ei=0, Le); sol1_030 := subs(e0=0, e1=1/3, ei=0, sol1); # change-of-variables transfo:=proc(L,a) local i,f; global x,Dx; f:=add(mult( subs(x=a,coeff(L,Dx,i)) , (1/diff(a,x) * Dx)$i ),i=0..degree(L,Dx)); sort(collect(f/lcoeff(f,Dx),Dx,factor),Dx) end; Lf := transfo(L030, f); sol_f := subs(x=f, sol1_030); # is a solution of Lf # Copy the file http://www.math.fsu.edu/~hoeij/files/equiv in this folder so that # the following read command succeeds: read "equiv"; # Compute the map from "solutions of Lf" to "solutions of P2". TR := equiv(Lf, P2); # Apply that map to sol_f to produce a solution of P2 sol_P2 := eval(diffop2de(TR,y(x)), y(x) = sol_f); sol_P2 := collect(%, indets(%,function), factor); # Check that this is indeed a solution: simplify( eval(diffop2de(P2, y(x)), y(x) = sol_P2) ); # OK, we solved the equation P2. # Next we'll see if we can shorten the expression. # First, lets interchange e1 and ei. Then we have (0,0,1/3). # interchanging 1,infinity is done by the mobius transformation x -> x/(x-1). # So we should update f as follows: f := factor( f/(f-1) ); # Lets compute the solution corresponding to (0,0,1/3). # After that, we'll fiddle a bit with (0,0,1/3), changing it by some integers (that have to add # up to an even integer) to see if we get a shorter solution. for es in [ [0,0,1/3], [-1,1,1/3]] do # Compute the hypergeometric equation with exponent differences es[1],es[2],es[3] # and a solution: Lnew := subs(e0=es[1], e1=es[2], ei=es[3], Le); sol1_new := subs(e0=es[1], e1=es[2], ei=es[3], sol1); # Now apply change of variables x -> f to: # Lnew (the hypergeometric solution with exponent-differences 1,1,1/3), and # sol1_new (one of its solutions) Lf := transfo(Lnew, f); sol_f := subs(x=f, sol1_new); # Compute the map from "solutions of Lf" to "solutions of P2". TR := equiv(Lf, P2); # Apply that map to sol_f to produce a solution of P2 sol_P2 := eval(diffop2de(TR,y(x)), y(x) = sol_f): lprint("this is the solution we get with e0,e1,e_infinity equal to", op(es)); sol_P2 := collect(%, indets(%,function), factor): od;