I'm trying to use Sage to find the general series solution to $y^{(4)}=\frac{y'y''}{1+x}$. So far my best efforts to derive the coefficient recurrence relations, inspired by a good book draft, have been along these lines:
R10 = QQ['a0,a1,a2,a3,a4,a5,a6,a7,a8,a9,a10']
a0,a1,a2,a3,a4,a5,a6,a7,a8,a9,a10 = R10.gens()
R.<x> = PowerSeriesRing(R10)
y = a0 + a1*x + a2*x^2 + a3*x^3 + a4*x^4 + a5*x^5 + \
a6*x^6 + a7*x^7 + a8*x^8 + a9*x^9 + a10*x^10 + O(x^11)
y1 = y.derivative()
y2 = y1.derivative()
y3 = y2.derivative()
y4 = y3.derivative()
f = (1+x)*y4-y1*y2
i = ideal(f)
# g = i.groebner_fan(); g.reduced_groebner_bases() # a wish
# q = R.quotient(i) # works, but not so useful by itself
and my other approaches ended in tracebacks:
x = var('x'); y = function('y', x)
desolve(diff(y,x,4)-diff(y,x)*diff(y,x,2)/(1 + x), y, contrib_ode=True)
# NotImplementedError: Maxima was unable to solve this ODE.
desolve_laplace(diff(y,x,4) - diff(y,x)*diff(y,x,2)/(1 + x), y)
# TypeError: unable to make sense of Maxima expression
I would like to at least solve that and determine the radius of convergence. Ideally (more generally), it would be nice to have a good bag of tricks for working with series DEs such as I imagine others have already created. For this, I would like to find or develop techniques to incorporate:
- more convenient coefficients, e.g. from this thread
- a way to derive the recurrence relations using Python's lambda operator or this nice trick
- a solver for higher order ODEs such as above
Any hints, links, references or suggestions would be appreciated.
