From: israel@math.ubc.ca (Robert Israel) Subject: Re: Nonlinear ODE Date: 19 Mar 2000 07:55:38 GMT Newsgroups: sci.math Summary: [missing] In article <8b1cn1$mr0$1@nnrp1.deja.com>, wrote: >Does anyone know of any nontrivial solutions to the following ODE > >y'' + y' + y^2 = 0 > >Does anyone know if this equation has a name? Maple's "odeadvisor" function calls this a modified Emden equation. And the solution is given as / d \ 2 y(x) = _a &where [{|--- _b(_a)| _b(_a) + _b(_a) + _a = 0}, \d_a / d {_a = y(x), _b(_a) = -- y(x)}, dx / | 1 {y(x) = _a, x = | ------ d_a + _C1}] | _b(_a) / Robert Israel israel@math.ubc.ca Department of Mathematics http://www.math.ubc.ca/~israel University of British Columbia Vancouver, BC, Canada V6T 1Z2 ============================================================================== From: rusin@vesuvius.math.niu.edu (Dave Rusin) Subject: Re: Why nobody answer this question? Date: 24 Oct 2000 08:05:32 GMT Newsgroups: sci.math I don't have a full answer to the original question, but I find this to be an interesting example, as I will explain below. Originally, in article , Zakir F. Seidov wrote: >Given ODE: > >(x^2 y')' = -x^2 y^2, > >with initial values: > >y(0)=1, and y'(0)=0. > >Series expansion (get by Mathematica) gives, for the convergence radius, >the numerical value 3.9626. What is the meaning of this value? In response, in article <8sqfit$qi1$1@nntp.itservices.ubc.ca>, Robert Israel wrote (correctly, of course): > >The solution in the complex plane has some sort of singularity >on the circle of radius 3.9626 centred at the origin, if that's the >radius. Actually the series, according to Maple, is > > 2 4 11 6 8 97 10 > y(x) = 1 - 1/6 x + 1/60 x - ---- x + 1/8505 x - -------- x + > 7560 10692000 [and more, deleted by djr] > >and assuming this pattern of signs continues indefinitely the >closest singularities to the origin will be on the imaginary axis. Taking a cue from this expansion we try replacing the (complex) variable x with u=-x^2; positive real values of u correspond to the points x on the imaginary axis. Using the chain rule we find the relationship between y and u to be / 2 \ /d \ |d | 2 6 |-- y(u)| + 4 u |--- y(u)| = y(u) \du / | 2 | \du / and the initial conditions y(0)=1, y'(0)=1/6. This ODE appears to admit the series solution y(u) = 1 + (1/6) u + (1/60) u^2 + ... consistent with Robert Israel's calculation. Furthermore, the ODE may be solved numerically without much trouble for small u (after using, say, the series solution to estimate y(eps) and y'(eps) for some small non-zero eps so as to avoid a division-by-zero error). For example, I find the solution curve has a point near (u, y) = (10, 10.2231422136984752) (where dy/du = 3.81661668823352369 or so). Now, it is proposed that there is a singularity of some kind near u = 3.9626^2 = 15.702, so we try to continue the numerical solution to larger u. This proves difficult, as indeed it appears there is a pole near this value of u, making both y and dy/du large. At this point it becomes preferable to study the curve on its side, that is, to think of u as a function of y, let y increase to infinity, and study the limiting value of u. In an attempt to study rates of convergence of u(y), I noticed my estimates were better expressed in terms of powers of z = y^(-1/2), which is more convenient anyway, since letting y approach infinity is equivalent to dropping z to the finite value 0. So let me do this _first_ (that is, before inverting the roles of u and the dependent variable). In the previous ODE substitute y = 1/z^2; then z = z(u) satisfies this equation: / 2 \ /d \ /d \2 |d | -12 z(u) |-- z(u)| + 24 u |-- z(u)| - 8 z(u) u |--- z(u)| = 1 \du / \du / | 2 | \du / Now the initial conditions are z(0)=1, z'(0)=-1/12. The previously found point on the curve now has coordinates (u, z) = (10, 0.312757547443019738) or thereabouts, and the tangent line has dz/du = -.0583810559410273223. We want to follow the curve to the point where z = 0 and expect to find there a point with u-coordinate near 15.7 . (One approach is to solve this ODE in series: z(u) = 1 - u/12 + ... and then solve the equation: [this approximation] = 0. I won't do it this way.) As I said before, it seems preferable to me to treat z as the independent variable, since the point we seek has a known z-coordinate (=0). To invert the equation, use the inverse function theorem: (dz/du) = 1 / (du/dz) and (d^2z/du^2) = - (d^2u/dz^2) / (du/dz)^3 . If I have done this correctly, the ODE is now: / 2 \ /d \2 /d \ |d | /d \3 -12 z |-- u(z)| + 24 u(z) |-- u(z)| + 8 z u(z) |--- u(z)| - |-- u(z)| \dz / \dz / | 2 | \dz / \dz / My numerical computations show that if we proceed from the point where (z, u) = (0.31... , 10) and du/dz = 1/(-.058...) , we terminate at z=0 with u = 15.7179392534512642 (and du/dz = -19.4224236922900625) (all approximate, of course). (Indeed, one may simply invert formally the previous series solution to get u = -12(z-1) + ... ; the coefficients of this series appear to decrease fairly rapidly so we may evaluate this series at (z-1) = -1 . Fifty terms gives the same value for u to within the magnitude of the last summand, about 10^(-12). ) Now, here is where I find the situation very interesting. The last ODE I wrote down now appears to have a solution which we find has u(0)=15.7... and u'(0)=-19.4... . These two numbers (nearly) satisfy the relation u(0) = u'(0)^2 / 24, which is convenient, since the ODE can be solved for u''(z) and we find a fraction with a numerical value at z=0, plus a term ( (u')^2 - 24 u ) / z which would mean u is not twice differentiable at z=0 unless u(0) = u'(0)^2 / 24 precisely. If we begin with that premise, we can attempt to find a series solution for the ODE near z = 0: for any nonzero constant c we find u(z) = 3 5 6 2 2 72 z 132192 z 12939264 z 7 1/24 c + c z + 18/5 z - -- ---- + ------ ---- - -------- ---- + O( z ) 25 c 625 3 3125 4 c c is a solution to the differential equation ... HOWEVER ... this series cannot be continued! Irrespective of any multiples of z^7, z^8, ... added to the terms shown, we find that when the series is substituted for u(z) in the ODE, there is a nonvanishing term beginning ( - 155271168/3125 / c^3 ) z^6 + ... I guess this shows there is no analytic solution u = u(z) to the differential equation near z=0, and further suppose this means that there is _not_ a simple pole but rather an essential singularity in the solution of the initial ODE, but I am rusty enough on my complex analysis that I would not put a lot of confidence in my conclusions here. Anyway, my computations lead me to believe the radius of convergence is closer to 3.9645856345211241 than to the 3.9626 mentioned by the original poster. I know of no significance to this number. Nothing nearby seemed appropriate in the tables at http://www.cecm.sfu.ca/projects/ISC/ISCmain.html I must say the original differential equation turned out to be more interesting than it first appeared. dave PS - This is almost, but not quite, an example of a Sturm-Liouville equation (p(x) y')' + q(x) y = 0. I know nothing about them, really, but perhaps some of the tools in S-L theory can be applied here too. ============================================================================== From: seidov@bgumail.bgu.ac.il (Zakir F. Seidov) Subject: Re: Why nobody answer this question? Date: 26 Oct 2000 15:04:28 -0400 Newsgroups: sci.math To Robert Israel, Dave Rusin , and all. ============= ========== === Now it's may be considered as established, thanks to your messages (and my own calcs), that for ODE LEE2: (x^2 y')'=-x^2 y^n, (n=2), y(0)=1,y'(0)=0, that series expansion of y(x) has radius of convergence (RC) (and Pade approximation of y(ix) has pole) at x=3.96458563... And I'm really very grateful for you for your REALLY VERY GREAT messages! Now the question is left (and this was the main drive of my first post on the LEE2): How is the radius of covergence of the serial expansion ============================== of the particular solution governed by the ODE - and in this particular case, =================== by value of n (taking initial values the same). In my old paper, in Sov.Astron.-AJ, vol. 21, pp 399-400 (1977) (yes! 23 years ago, and note that this is Engl. transl. of the more early original paper in Russian), I wrote: "It is intuitevely clear that the radius of convergence ... of the series... is a decreasing function of the index n... Unfortunately, no simple means exists for analyzing the convergence of the series, the coefficients of which are given in a rather cumbersome form by the recursion relation (12)" The recursion relation for coefficients of the serial expansion of the particular solution of the ODE LEEn (as I call it) is: y(x)=1 +a_1 x^2+...+a_i x^(2k)+... a^(m+1)=1/((m)(m+1)(2m+3))* Sum_{i=1}^{i=m}[(i*n+i-m)(m-i+1)(3+2(m-i))a_i*a_(m+1-i)], m>0, a_1=-1/6, a_2=n/120, a_3=-n(8n-5)/3*7!, a_4=n(122n^2-183n+70)/9*9!, and so on. This relation does greatly reduce the task of calculation of the series' coefficients, of RC, of Pade approximation, of "pole' of y(ix). And besides, I hope, it may help to some expert matematician to find RC analytically from recurrence relation for coefficicents. ============ Note, e.g. that as it may be shown, for integer n, the series has coefficients of alternate signs. For general background: ======================= This equation describe the internal structure of polytropic stars, and physically reasonable values of n are in the interval [o,5]. Also LEEn (with the same initial values!) has known solutions: at n=0, y(x)=1-1/6 x^2 (and recurrence relation confirms this!), at n=0, y(x)=sin(x)/x (and recurrence relation confirms this!, and in this case RC=Infinity ), at n=5, y(x)=(1+1/3 x^2))^(-1/2), (and recurrence relation confirms this!, and in this case RC=3^(1/2)). ZFS ============================================================================== From: israel@math.ubc.ca (Robert Israel) Subject: Re: Why nobody answer this question? Date: 24 Oct 2000 19:13:26 GMT Newsgroups: sci.math In article <8sqhcl$9s4@mcmail.cis.McMaster.CA>, Zdislav V. Kovarik wrote: >In article <8sqfit$qi1$1@nntp.itservices.ubc.ca>, >Robert Israel wrote: >:In article , >:Zakir F. Seidov wrote: >:>Given ODE: >:>(x^2 y')' = -x^2 y^2, >:>with initial values: >:>y(0)=1, and y'(0)=0. >:>Series expansion (get by Mathematica) gives, >:>for the convergence radius, >:> the numerical value 3.9626. >:>What is the meaning of this value? >:The solution in the complex plane has some sort of singularity >:on the circle of radius 3.9626 centred at the origin, if that's the >:radius. Actually the series, according to Maple, is >: 2 4 11 6 8 97 10 >: y(x) = 1 - 1/6 x + 1/60 x - ---- x + 1/8505 x - -------- x + >: 7560 10692000 >: 457 12 98239 14 13339 16 >: --------- x - ------------- x + ------------- x - >: 673596000 1980372240000 3740703120000 >: 3276433 18 20 >: ----------------- x + O(x ) >: 12953119728780000 >:and assuming this pattern of signs continues indefinitely the >:closest singularities to the origin will be on the imaginary axis. >Out of curiosity: will a Pade expansion suggest such a singularity? Taking the series to O(x^50) and doing a Pade expansion of degree (24,24), Maple says the denominator has roots at 3.9645863252023634436 I and 3.9645849438502091794 I (and their complex conjugates). Presumably this actually indicates second-order poles of the solution. In fact, it's not hard to see that any pole of a solution of the differential equation must be second-order, and have principal part -6/(x-a)^2 + 12/(5 a (x-a)) (if the pole is at x=a). However, I don't know whether the solutions are meromorphic in the complex plane. Robert Israel israel@math.ubc.ca Department of Mathematics http://www.math.ubc.ca/~israel University of British Columbia Vancouver, BC, Canada V6T 1Z2 ============================================================================== From: israel@math.ubc.ca (Robert Israel) Subject: Re: Why nobody answer this question? Date: 25 Oct 2000 18:07:07 GMT Newsgroups: sci.math In article <8t4n0m$fud$1@nntp.itservices.ubc.ca>, Robert Israel wrote: >Taking the series to O(x^50) and doing a Pade expansion of degree (24,24), >Maple says the denominator has roots at 3.9645863252023634436 I and >3.9645849438502091794 I (and their complex conjugates). Presumably this >actually indicates second-order poles of the solution. In fact, it's not >hard to see that any pole of a solution of the differential equation must >be second-order, and have principal part >-6/(x-a)^2 + 12/(5 a (x-a)) >(if the pole is at x=a). Of course, I was assuming a <> 0 here. There is a solution with a pole at x=0, namely y = -2/x^2. > However, I don't know whether the solutions are >meromorphic in the complex plane. I agree with Dave Rusin: there isn't really a pole at x=a. The expansion continues y = -6/(x-a)^2 + 12/(5 a (x-a)) - 48/(25 a^2) + 204/(125 a^3) (x-a) - 876/(625 a^4) (x-a)^2 + 3612/(3125 a^5) (x-a)^3 but then you need a logarithmic term to keep going. One of many interesting transformations you can do to this Emden equation is x = exp(-t), y = exp(2 t) u, which makes it into an overdamped nonlinear spring: u'' + 3 u' + 2 u + u^2 = 0 The equilibrium solutions u=0 and u=-2 correspond to y=0 and y = -2/x^2. The equilibrium u=0 has eigenvalues -1 and -2. Zakir's series y = 1 - x^2/6 + ... corresponds to a solution approaching the equilibrium u=0 in the phase plane along the eigenvector for eigenvalue -2 with u > 0. The solution approaching in the opposite direction is y = -1 - x^2/6 - ... (if Zakir's solution is y1(x), this one is -y1(i x)). The radius of convergence of the power series corresponds to this solution going off to infinity at a finite value of t. A look at the phase plane seems to indicate that except for the solutions coming out of the saddle point u=-2 as t -> -infinity, all the rest will go off to infinity in some direction as t decreases (probably to a finite value rather than -infinity). Robert Israel israel@math.ubc.ca Department of Mathematics http://www.math.ubc.ca/~israel University of British Columbia Vancouver, BC, Canada V6T 1Z2