Ask Your Question
2

series solutions of higher order ODEs

asked 2012-02-15 01:33:26 -0500

bgins gravatar image

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.

edit retag flag offensive close merge delete

1 answer

Sort by ยป oldest newest most voted
0

answered 2017-05-26 13:20:22 -0500

dan_fulea gravatar image

The substitution $z=y'$ is natural, so let us please then forget about the posted $y$.

We will not use $z$. Let us use $y$ "back" instead of $z$ in the sequel. So we solve (instead) around $x=0$: $$ y'' ' = y\ y' \ /\ (1+x)\ .$$

Reusing the above code...

R10.<a0,a1,a2,a3,a4,a5,a6,a7,a8,a9,a10> = \
    QQ['a0,a1,a2,a3,a4,a5,a6,a7,a8,a9,a10']

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()

f = (1+x)*y3 - y*y1

# J = ideal(f) # --> leads to an error
f.list()

... we get the coefficients of $f$:

[-a0*a1 + 6*a3,
 -a1^2 - 2*a0*a2 + 6*a3 + 24*a4,
 -3*a1*a2 - 3*a0*a3 + 24*a4 + 60*a5,
 -2*a2^2 - 4*a1*a3 - 4*a0*a4 + 60*a5 + 120*a6,
 -5*a2*a3 - 5*a1*a4 - 5*a0*a5 + 120*a6 + 210*a7,
 -3*a3^2 - 6*a2*a4 - 6*a1*a5 - 6*a0*a6 + 210*a7 + 336*a8,
 -7*a3*a4 - 7*a2*a5 - 7*a1*a6 - 7*a0*a7 + 336*a8 + 504*a9,
 -4*a4^2 - 8*a3*a5 - 8*a2*a6 - 8*a1*a7 - 8*a0*a8 + 504*a9 + 720*a10]

We give or consider a0, a1, a2 as given.

Solving by the Cauchy-Kovalevskaya "procedure", we have to determine recursively a3, a4, a5, a6, a7, a8, a9, a10, and so on, starting from a0, a1, a2, by making each line evaluate to $0$. We get manually $$ a_3 = \frac 16 a_0a_1\ , $$ $$ a_4 = \frac 1{24}( a_1^2 + 2a_0a_2 + a_0a_1 )\ , $$ and so on.

I think it is not a good idea to isolate in this way all these equations, then start solving. Instead, let us write the given DE in the form $$ y'' ' = y\ y'\ (1-x+x^2-x^3+\dots)\ , $$ and consider it formally in the ring $\mathbb Q[[x]]$. Then the equations in $a_3,a_4,\dots$ come already separated, and we can simply write the code to get the first some $100$ coefficients.

(I am not willing to write optimal code, time...)

To have a precise situation, let us better use values instead of the symbols $a_0,a_1, a_2$. For instance $a_0=0$, $a_1=1$, $a_2=2$. Then...

a = [0,1,2]    # use other initial conditions, next time
PSR.<x> = PowerSeriesRing( QQ, 100 )
C = 1 / PSR(1+x)

for k in [ 3..99 ]:
    print k,
    # trying to get a(k) from y = sum( [ a(k) * x^k for k in [0..ooh] ] )
    yk = sum( [ a[j] * x^j for j in range(k) ] ) + O(x^k)

    ak = 1/(k*(k-1)*(k-2)) * ( yk * diff( yk, x ) * ( C+O(x^k) ) ).list()[k-3]
    print ak
    a.append( ak )

y = PSR( sum( [ a[k] * x^k for k in [0..99] ] ) )

This computes some values. For instance, the first some $30$ entries are:

sage: a[:30]
[0,
 1,
 2,
 0,
 1/24,
 1/12,
 1/40,
 -67/5040,
 13/1152,
 -7/1440,
 4621/1209600,
 -38399/13305600,
 11287/4838400,
 -12737/6988800,
 2351057/1614412800,
 -171374993/145297152000,
 1692465871/1743565824000,
 -115101113/142502976000,
 32202455561/47424990412800,
 -4213094929/7313918976000,
 22201339257697/45053740892160000,
 -267905698350809/630752372490240000,
 5115211318062343/13876552194785280000,
 -8562348061276901/26596725040005120000,
 928302494123858617/3282795776366346240000,
 -7613558036179561/30490363247984640000,
 1985785680757233642821/8962032469480125235200000,
 -590031815878576041977/2987344156493375078400000,
 33273819697787429364371/188202681859082629939200000,
 -3468510498746975345435917/21831511095653585072947200000]

(These are exactly computed.) Let us check the solution modulo $O(x^{100})$:

# : diff( y,x,3 ) - y * diff(y,x) / (1+x)
# the above is a mess with terms in degrees 97, 98, 99, but..
sage: diff( y,x,3 ) - y * diff(y,x) / (1+x) + O( x^97 )
O(x^97)

The radius of convergence is in this case probably one, one can try...

for k in [1..99]:
    ak = a[k]
    print ( "a[%2d] ~ %13.10f and a[%2d]^(1/%2d) ~ %8.6f"
            % ( k, RR(ak), k, k, abs(RR(ak))^(1/k) ) )

and the first and last lines shown are:

a[ 1] ~  1.0000000000 and a[ 1]^(1/ 1) ~ 1.000000
a[ 2] ~  2.0000000000 and a[ 2]^(1/ 2) ~ 1.414214
a[ 3] ~  0.0000000000 and a[ 3]^(1/ 3) ~ 0.000000
a[ 4] ~  0.0416666667 and a[ 4]^(1/ 4) ~ 0.451801
a[ 5] ~  0.0833333333 and a[ 5]^(1/ 5) ~ 0.608364
a[ 6] ~  0.0250000000 and a[ 6]^(1/ 6) ~ 0.540742
a[ 7] ~ -0.0132936508 and a[ 7]^(1/ 7) ~ 0.539447
a[ 8] ~  0.0112847222 and a[ 8]^(1/ 8) ~ 0.570902
a[ 9] ~ -0.0048611111 and a[ 9]^(1/ 9) ~ 0.553313
:::::::::::::::::::::::::::::::::::::::::::::::::
a[90] ~  0.0000051223 and a[90]^(1/90) ~ 0.873406
a[91] ~ -0.0000049542 and a[91]^(1/91) ~ 0.874386
a[92] ~  0.0000047934 and a[92]^(1/92) ~ 0.875348
a[93] ~ -0.0000046394 and a[93]^(1/93) ~ 0.876295
a[94] ~  0.0000044920 and a[94]^(1/94) ~ 0.877225
a[95] ~ -0.0000043507 and a[95]^(1/95) ~ 0.878140
a[96] ~  0.0000042154 and a[96]^(1/96) ~ 0.879040
a[97] ~ -0.0000040855 and a[97]^(1/97) ~ 0.879925
a[98] ~  0.0000039610 and a[98]^(1/98) ~ 0.880796
a[99] ~ -0.0000038414 and a[99]^(1/99) ~ 0.881653

I have maybe to stop here with sage sample code.

About the convergence radius... The recurrence relation is thus of the shape $a_k=$ a sum of monomials of order two in the previous variables, with coefficients that can be recursively "controlled". With some work, one may (??) maybe show that the radius of convergence is one. (The one is my quick guess, based on the fact that the limit $r=\lim |a_k|^{1/k}$ may / should exist, and then we expect $r\sim r^2$. Minoration / majoration for the "polynomial part in $k$" is constricted by the power $1/k$. But it is just a guess.)

Note: The use of solve_laplace really wants to perform a Laplace transform. But which is the first step when we try to apply it on the right side of the given DE ?!

Note: The answer is not given mot-a-mot to the posted (relatively general) question, but in the sense of the progress. (Which is subjective.) I think, sample code serves good as an answer.

edit flag offensive delete link more

Your Answer

Please start posting anonymously - your entry will be published after you log in or create a new account.

Add Answer

Question Tools

3 followers

Stats

Asked: 2012-02-15 01:33:26 -0500

Seen: 332 times

Last updated: May 26