ASKSAGE: Sage Q&A Forum - RSS feedhttps://ask.sagemath.org/questions/Q&A Forum for SageenCopyright Sage, 2010. Some rights reserved under creative commons license.Fri, 26 May 2017 20:20:22 +0200series solutions of higher order ODEshttps://ask.sagemath.org/question/8717/series-solutions-of-higher-order-odes/I'm trying to use Sage to find the general series solution to [$y^{(4)}=\frac{y'y''}{1+x}$](http://math.stackexchange.com/questions/109189/basic-reference-material-about-odes-such-as-saparability-with-calculations-and-e "a posted math question"). So far my best efforts to derive the coefficient recurrence relations, inspired by a good [book draft](http://wdjoyner.com/teach/DiffyQ/des-book.pdf "David Joyner (2007), Introductory Differential Equations using SAGE"), 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](http://ask.sagemath.org/question/441/writing-re-usable-sage-scripts#785 "Writing re-usable sage scripts") have [already](file:///home/mike/Projects/Sage/Differential_Equations/series.sage "Are you the right mike?") created. For this, I would like to find or develop techniques to incorporate:
- more convenient coefficients, e.g. from [this thread](http://ask.sagemath.org/question/611/implicitly-defining-a-sequence-of-variables "implicitly-defining-a-sequence-of-variables")
- a way to derive the recurrence relations using Python's lambda operator or this [nice trick](http://ask.sagemath.org/question/142/generating-series#308 "LazyPowerSeriesRing.product_generator.compute_coefficients")
- a solver for higher order ODEs such as above
Any hints, links, references or suggestions would be appreciated.Wed, 15 Feb 2012 08:33:26 +0100https://ask.sagemath.org/question/8717/series-solutions-of-higher-order-odes/Answer by dan_fulea for <p>I'm trying to use Sage to find the general series solution to <a href="http://math.stackexchange.com/questions/109189/basic-reference-material-about-odes-such-as-saparability-with-calculations-and-e" title="a posted math question">$y^{(4)}=\frac{y'y''}{1+x}$</a>. So far my best efforts to derive the coefficient recurrence relations, inspired by a good <a href="http://wdjoyner.com/teach/DiffyQ/des-book.pdf" title="David Joyner (2007), Introductory Differential Equations using SAGE">book draft</a>, have been along these lines:</p>
<pre><code>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
</code></pre>
<p>and my other approaches ended in tracebacks:</p>
<pre><code>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
</code></pre>
<p>I would like to at least solve that and determine the radius of convergence. <em>Ideally</em> (more generally), it would be nice to have a good bag of tricks for working with series DEs such as I imagine <a href="http://ask.sagemath.org/question/441/writing-re-usable-sage-scripts#785" title="Writing re-usable sage scripts">others</a> have <a title="Are you the right mike?">already</a> created. For this, I would like to find or develop techniques to incorporate:</p>
<ul>
<li>more convenient coefficients, e.g. from <a href="http://ask.sagemath.org/question/611/implicitly-defining-a-sequence-of-variables" title="implicitly-defining-a-sequence-of-variables">this thread</a></li>
<li>a way to derive the recurrence relations using Python's lambda operator or this <a href="http://ask.sagemath.org/question/142/generating-series#308" title="LazyPowerSeriesRing.product_generator.compute_coefficients">nice trick</a></li>
<li>a solver for higher order ODEs such as above</li>
</ul>
<p>Any hints, links, references or suggestions would be appreciated.</p>
https://ask.sagemath.org/question/8717/series-solutions-of-higher-order-odes/?answer=37709#post-id-37709The 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.
Fri, 26 May 2017 20:20:22 +0200https://ask.sagemath.org/question/8717/series-solutions-of-higher-order-odes/?answer=37709#post-id-37709