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;