read("hypergeomdeg3");
# The differentiation variable is hard-coded as x in hypergeomdeg3
L2r := subs({t=x,Dt=Dx,x=y},L2):
# Hypergeomdeg3 is not made to work efficiently with rational fraction
# coefficients, so instead of forcing it, we use evaluation and
# interpolation to reduce the problem to that case.
### Evaluation: substitute different values for y in the differential
### operator and find a hypergeometric solution with a degree 3
### substitution using hypergeomdeg3
nvals := 15:
pts := []:
vals := []:
for yy from 11 to 10+nvals do
L2_y := subs(y=yy,L2r):
val := simplify(hypergeomdeg3(L2_y,x)):
vals := [op(vals),val[1]]:
pts := [op(pts),yy]:
od:
### Interpolation and rational reconstruction:
### We want to find fractions in y for each coefficient in the
### hypergeometric expressions we have found. So first we need to
### separate those expressions into coefficients, or, rather, convert
### them into polynomials and let Maple identify their coefficients.
# First we identify the hypergeoms and the radicals which cannot be
# polynomial variables
HR := []:
for i from 1 to nvals do
hr := indets(vals[i],{function,radical}):
HR := [op(HR),hr]:
od:
# Verify that they are not changing for different points
print(HR):
# If necessary, throw away the points for which they are different
# vals := subsop(1=NULL, vals):
# pts := subsop(1=NULL, pts):
HR := HR[1]:
# So we have one radical and two hypergeoms. We can clear out the
# radical by multiplying and replace the hypergeoms with polynomial
# variables.
# Note that the two hypergeoms are linearly independent and only
# appear with degree 1 in the expressions, so it is acceptable to
# treat them as polynomial variables.
# We collect what remains in lists of polynomials (numerator and
# denominator) for interpolation
nums := []:
dens := []:
for i from 1 to nops(pts) do
tmp := 1/HR[1] * vals[i]:
tmp := algsubs(HR[2]=H1, tmp):
tmp := algsubs(HR[3]=H2, tmp):
num := numer(tmp):
num := 1/lcoeff(num) * num:
den := denom(tmp):
den := 1/lcoeff(den) * den:
nums := [op(nums),num]:
dens := [op(dens),den]:
od:
# Now we can do the interpolation
num_rec1 := CurveFitting[PolynomialInterpolation](pts,nums,y):
den_rec1 := CurveFitting[PolynomialInterpolation](pts,dens,y):
# And reconstruct a rational fraction from the result
m := mul(y-yy, yy in pts):
num_rec := ratrecon(num_rec1,m,y):
den_rec := ratrecon(den_rec1,m,y):
SOL := num_rec/den_rec;
# Now we substitute back and try to simplify the expression.
# Simplify the hypergeoms
hyper1 := (1-27*x^3)^(-1/3) * hypergeom([-1/3,-2/3],[1],27*x^3); # Pfaff
hyper2 := (1-27*x^3)^(2/3) * hypergeom([-2/3,2/3],[2],27*x^3); # Pfaff
# Check
series(HR[2] - hyper1, x=0, 50);
# 0
series(HR[3] - hyper2, x=0, 50);
# 0
SOL := subs({H1=hyper1,H2=hyper2}, SOL);
SOL := HR[1]*SOL;
SOL := subs({x=t,y=x},SOL);
normal(SOL);
# Observation
SOL := map(collect,factor(SOL),indets(SOL,{function}),factor);
# There are some cubic roots of 27^t^3-1 which maple can't simplify.
# But we know that they are equal up to a constant (cubic root of 1),
# and since we are only finding the solution up to constant, we might
# as well choose the factor to be -1.
inds := indets(SOL,{radical});
SOL := (1-27*t^3)^(-1/3) * 1/inds[1] * inds[2] * SOL;
# Check
simplify(add(coeff(L2,Dt,i)*diff(SOL,t$i),i=1..2)+coeff(L2,Dt,0)*SOL);
# 0
# Remove the leading coefficient
SOL := 81*SOL:
lprint(map(collect,factor(SOL),indets(SOL,{function}),factor));
# Here the obstinate cubic root is gone. Let's just keep this
# expression instead of keeping the unsimplified one.
SOL_old := SOL:
SOL := -((6*t^2*x^3+3*t^2+3*t*x-2*x^2)*hypergeom([-2/3, -1/3],[1],27*t^3)+5*(2*t*x^3+t-x)*t*hypergeom([-2/3, 2/3],[2],27*t^3))*(x^3+2)*(4*x^3-1)/(2*x^3+1)/t^3/(t*x^3+2*t-x)/(4*t^2*x^3-t^2+2*t*x-x^2): # printed above
# Maple can verify that the two solutions are equal
simplify(SOL-SOL_old);
# 0
# Remove the factors constant in t
SOL := SOL/((x^3+2)*(4*x^3-1)/(2*x^3+1));
# Final check
simplify(add(coeff(L2,Dt,i)*diff(SOL,t$i),i=1..2)+coeff(L2,Dt,0)*SOL);
# 0
# Local Variables:
# mode: maplev
# End: