From: rusin@vesuvius.math.niu.edu (Dave Rusin) Subject: Re: Differential Equation Problem Date: 17 Nov 2001 00:48:49 GMT Newsgroups: sci.math Summary: Sample algebraic approaches for solving an autonomous system of ODEs In article <3BF54EFC.B8E7ADE1@NOSPAMhobonet.com>, Jake wrote: >dy/dx = C1 - K*y*Sqrt[y^2 + z^2] > >dz/dx = C2 - K*z*Sqrt[y^2 + z^2] > >Would someone please tell me if the above equations have an exact >solution or must I use numerical methods? Well, I don't know about "must". You can try to make progress toward a symbolic solution. Take advantage of the "circle-like" form of the equations -- here's a suggestion in Maple-ese: E1:=diff(y(x),x) = C1 - K*y(x)*sqrt(y(x)^2 + z(x)^2); E2:=diff(z(x),x) = C2 - K*z(x)*sqrt(y(x)^2 + z(x)^2); subs({y(x)=r(x)*(1-m(x)^2)/(1+m(x)^2), z(x)=r(x)*2*m(x)/(1+m(x)^2)},[E1,E2]): simplify(%,symbolic); aha:=solve({op(%)},{diff(r(x),x),diff(m(x),x)}); You'll get a new pair of differential equations of the form dr/dx = ..., dm/dx = ... One of them is this: 2 d C2 m(x) + 2 C1 m(x) - C2 -- m(x) = - 1/2 -------------------------, dx r(x) which is readily solved for r, so we can eliminate this variable altogether: solve(aha[1],r(x)); subs(r(x)=%,aha[2]); factor(lhs(%)-rhs(%)); de:=numer(%); The result is a _single_ differential equation in _one_ unknown: /d \2 /d \2 3 /d \2 0 = -12 C2 m(x) |-- m(x)| - 4 |-- m(x)| C2 m(x) - 8 C1 |-- m(x)| \dx / \dx / \dx / / 2 \ / 2 \ / 2 \ |d | 4 |d | |d | 3 + 2 |--- m(x)| C2 m(x) + 4 |--- m(x)| C1 m(x) + 4 |--- m(x)| C1 m(x) | 2 | | 2 | | 2 | \dx / \dx / \dx / / 2 \ |d | 6 2 5 4 2 - 2 |--- m(x)| C2 + K m(x) C2 + 4 K m(x) C2 C1 - K m(x) C2 | 2 | \dx / 4 2 2 2 2 2 2 + 4 K m(x) C1 - K m(x) C2 + 4 K C1 m(x) - 4 K C1 m(x) C2 + K C2 Well, is this one solvable? I threw it to Maple, and it churned out a big mess. From that mess I discern that it is better to solve for _x_ in terms of _m_ and then invert. One can convert the equation as follows: subs(diff(m(x),x$2)=-xmm/xm^2, de); subs(diff(m(x),x)=1/xm, %); subs(m(x)=m, %); ans:=numer(-%)/(1+m^2); The result is interesting. Since x itself does not appear here, this equation can be expressed in the form 2 Quad w' - K*(Quad w)^2 + 4 Cub/(1+m^2) = 0 where w = dx/dm (and w' = d^2x/dm^2 ) and where 3 2 Cub = 3 C2 m + C2 m + 2 C1, Quad = C2 m + 2 C1 m - C2 Cub:=(C2*m^3+3*C2*m+2*C1): Quad:=(C2*m^2+2*C1*m-C2): Lin:=C2*m +C1: ans2:= 2*Quad*xmm - K*(Quad*xm)^2 + 4*Cub/(1+m^2); So this brings us down to a _single_ _first_-order ODE involving only rational functions of m. Presumably there is a Risch-like procedure to decide whether or not the solutions can be expressed in "finite terms", but I don't know any details. (I tried recasting the ODE in terms of v = (K/2) Quad w, so that it has the form ( v' - v^2 ) + (junk) v + (junk) = 0, but that didn't seem like enough of an improvement.) I'm afraid I didn't get any further with this equation, but it does sort of look like something the pros can do. >What is your best method of solution? "Best"? Your choice. If you need numbers, solve numerically. dave ============================================================================== From: rusin@vesuvius.math.niu.edu (Dave Rusin) Subject: Re: Differential Equation Problem Date: 17 Nov 2001 19:15:59 GMT Newsgroups: sci.math In article <3BF54EFC.B8E7ADE1@NOSPAMhobonet.com>, Jake wrote: >dy/dx = C1 - K*y*Sqrt[y^2 + z^2] > >dz/dx = C2 - K*z*Sqrt[y^2 + z^2] > >Would someone please tell me if the above equations have an exact solution Here's a different way to think about that pair of equations. Since x never occurs explicitly in your system, this is an "autonomous" system of differential equations. Think of x as representing time, and y and z as variables which change in time. In my previous post I looked to see whether y and z could be given as some explicit function of t (it looks like the answer is "no"). But maybe you don't need that information. Maybe it's enough to know what combinations of values for y and z can ever occur together, without really paying attention to _when_ a particular combination can occur. That is, if, for example, it should happen to pass that the solution to a similar system of ODEs is, say, y = A sin( sqrt(x)+27/x ) z = A cos( sqrt(x)+27/x ) then you can spot right off the bat that the values of y and z always satisfy y^2 + z^2 = A^2 ; the solutions lie on circles at the origin. (Just which circle they fall on will depend on the initial conditions for the system.) How exactly you move around the circle in time is messy, but may be irrelevant. If that's what you want to know, then you can get a sense of what "orbits" the points (y,z) will traverse in your autonomous system as follows: if ever the orbit passes through a particular point (y,z) in the y-z plane, it does so with a slope equal to (dz/dx) / (dy/dx). That means the orbit will be a curve in the plane which looks like a the graph of a function z = z(y) which is a solution to the single ODE dz C2 - K*z*Sqrt[y^2 + z^2] (*) -- = ------------------------ dy C1 - K*y*Sqrt[y^2 + z^2] As a practical matter you could, for any given combination of C1, C2, K compute the right-hand side for a wide range of pairs (y,z); at each such point you can draw a short line segment with that slope to represent the tangent line of a solution curve passing through that point. Then let your eye join together the little marks into smooth curves. Those curves are the graphs of the functions z = z(y) which solve this ODE. As further proof of the proposition that algebraists should never be permitted to get anywhere near a differential equation, I will now consider the possibility of describing these curves in "closed form" for your particular system. That is, I will "solve" (*). First, I want to avoid those square roots, so I will re-write your initial system a little bit. In Maple, the steps I take are: E1:=diff(y(x),x) = C1 - K*y(x)*sqrt(y(x)^2 + z(x)^2); E2:=diff(z(x),x) = C2 - K*z(x)*sqrt(y(x)^2 + z(x)^2); simplify(subs({y(x)=u(x)^2-v(x)^2,z(x)=2*u(x)*v(x)},[E1,E2]),symbolic); solve({op(%)},{diff(u(x),x),diff(v(x),x)}); which gives us this square-root-free system to solve: 5 4 3 2 du -u C1 + K u + K u v - C2 v + 2 K u v -- = - 1/2 ---------------------------------------- dx 2 2 u + v 2 3 4 5 dv 2 K u v - u C2 + K u v + C1 v + K v -- = - 1/2 --------------------------------------- dx 2 2 u + v Then, as outlined above, I divide through to get a single ODE linking u and v: du/dv = (-u*C1+K*u^5+K*u*v^4-C2*v+2*K*u^3*v^2)/(2*K*u^2*v^3-u*C2+K*u^4*v+C1*v+K*v^5) As it happens, this ODE _does_ admit a "closed-form" solution, which is however much too painful to display. Maple must've gone through a number of substitutions to get the answer, and I'm probably undoing them to present the answer in better form. It looks better if we introduce a quantity R = sqrt( C1^2 + C2^2 ), and replace one of the variables (say v) with w = (C2*v/u + C1)/R, that is, we let v = u*(R*w-C1)/C2. In that case, the solution may be expressed as the curves of the form arctanh(w) = P1 + P2/u^4 + C3 where P1 and P2 are rational functions of w, and C3 is an arbitrary constant (used to distinguish the different solution curves). Precisely, it's P1:=w*(-2*R^2*w^2+w^2*C2^2+4*C1*R*w-2*R^2+C2^2)/(C2^2*(w^2-1)^2) P2:=C2^2/(2*R*K*(w^2-1)^2) So we can solve this algebraic equation for u: u = (P2/(arctanh(w)-P1-C3))^(1/4) If you prefer, you may use this to parameterize the curves in the u-v plane, and then in the y-z plane, using this w as a parameter, e.g. u = (P2/(arctanh(w)-P1-C3))^(1/4) v = u*(R*w-C1)/C2 Of course, you might at that point prefer a different parameter, e.g. use t=arctanh(w), so w=tanh(t)=(exp(t)-exp(-t))/(exp(t)+exp(-t)). Since the parameterized curves in the u-v or y-z plane will then involve both t and exp(t), it seems unlikely that these curves can be described in a pleasant implicit form without using, say, the Lambert-W function. On the other hand, you could take this parameterization and use it to rewrite the original ODE as a single first-order equation describing this parameter t as a function of x. In other words, you would not only know which points (y,z) you traverse, but you will know when (i.e. for which x) you get to each of them. I'll leave all these options to the reader. dave ============================================================================== From: "Peter L. Montgomery" Subject: Re: Differential Equation Problem Date: Sat, 17 Nov 2001 21:49:36 GMT Newsgroups: sci.math In article <9t6d1f$ior$1@news.math.niu.edu> rusin@vesuvius.math.niu.edu (Dave Rusin) writes: >In article <3BF54EFC.B8E7ADE1@NOSPAMhobonet.com>, >Jake wrote: > >>dy/dx = C1 - K*y*Sqrt[y^2 + z^2] >> >>dz/dx = C2 - K*z*Sqrt[y^2 + z^2] >> >>Would someone please tell me if the above equations have an exact solution > > >Here's a different way to think about that pair of equations. > >Since x never occurs explicitly in your system, this is an "autonomous" >system of differential equations. Think of x as representing time, >and y and z as variables which change in time. (text deleted) >As further proof of the proposition that algebraists should never >be permitted to get anywhere near a differential equation, I will >now consider the possibility of describing these curves in "closed form" >for your particular system. That is, I will "solve" (*). This algebraist (and number theorist) used maple too, and found that con1 is constant (i.e., independent of x) where |\^/| Maple V Release 5 (WMI Campus Wide License) ._|\| |/|_. Copyright (c) 1981-1997 by Waterloo Maple Inc. All rights \ MAPLE / reserved. Maple and Maple V are registered trademarks of <____ ____> Waterloo Maple Inc. | Type ? for help. > Cr := sqrt(C1^2 + C2^2); 2 2 1/2 Cr := (C1 + C2 ) > r := sqrt(y^2 + z^2); 2 2 1/2 r := (y + z ) > > con1 := Cr * (C1*(C1 - K*y*r) + C2*(C2 - K*z*r)) / (C2*y - C1*z)^2 > - K * arctanh((C1*y + C2*z)/(r*Cr)); con1 := 2 2 1/2 2 2 1/2 2 2 1/2 (C1 + C2 ) (C1 (C1 - K y (y + z ) ) + C2 (C2 - K z (y + z ) )) ------------------------------------------------------------------------ 2 (C2 y - C1 z) C1 y + C2 z - K arctanh(---------------------------) 2 2 1/2 2 2 1/2 (y + z ) (C1 + C2 ) > # Can replace arctanh by +- arcsinh((C1*y + C2*z)/(C2*y - C1*z)). > > # Verify that con1 is constant > > y := yf(x); y := yf(x) > z := zf(x); z := zf(x) > derivs := {diff(y, x) = C1 - K*y*r, diff(z, x) = C2 - K*z*r}; d 2 2 1/2 derivs := {-- zf(x) = C2 - K zf(x) (yf(x) + zf(x) ) , dx d 2 2 1/2 -- yf(x) = C1 - K yf(x) (yf(x) + zf(x) ) } dx > > factor(simplify(subs(derivs, diff(con1, x)))); bytes used=1015332, alloc=851812, time=0.25 0 > # If you investigate further, try using C1*y + C2*z and C2*y - C1*z # as your major variables, instead of y and z. > > ;quit; bytes used=1231164, alloc=851812, time=0.30 -- George W. Bush tells Americans not to be afraid to fly. Americans are flying their flags. Peter-Lawrence.Montgomery@cwi.nl Home: San Rafael, California Microsoft Research and CWI ============================================================================== From: rusin@vesuvius.math.niu.edu (Dave Rusin) Subject: Re: Differential Equation Problem Date: 19 Nov 2001 09:06:49 GMT Newsgroups: sci.math Jake wrote: > dy/dx = C1 - K*y*Sqrt[y^2 + z^2] > dz/dx = C2 - K*z*Sqrt[y^2 + z^2] ayatollah potassium wrote: > Let Z = y + i*z (complex number), C = C1+ i*C2, k = K > dZ/dx = C - k Z * |Z| (C complex, k real, x real) > > Rescaling x to Ax multiplies (C,k) by A (A real). > Rescaling Z to Bz changes (C,k) to (C/B, k*|B|), (B complex). ... and some other stuff which didn't quite work. Then, in article <3BF83D28.4C60174D@bibimus.edu>, s/he wrote: > I don't know if this substantially simplifies the computations > given in the other postings. It does. Let me take this opportunity to combine this with both my previous posts into a better framework. The rescalings allow us to assume the problem to solve is just dy/dx = 1 - y*Sqrt[y^2 + z^2] dz/dx = - z*Sqrt[y^2 + z^2] So it's really a _single_ system of DEs, not a family of them. As far as the phase-plane analysis goes (i.e. avoiding references to x), I still think this is most easily diagnosed with the substitutions y=u^2-v^2, z=2uv; then these new variables are linked by the single ODE dv/du = (v/u)*( (v+u)^2 + 1 )/( (v+u)^2 - 1 ) which is explicitly solveable; the solutions are the curves of the form 4 ln(v/u) + (v/u)^2 - (u/v)^2 + 1/(uv)^2 = C and, if necessary, this may be expressed in terms of y,z . These curves are easily parameterized using w = v/u : u = w^(-1/2)/( C - 4*ln(w) - w^2 + 1/w^2 )^(1/4) v = w^(+1/2)/( C - 4*ln(w) - w^2 + 1/w^2 )^(1/4) This lets us parameterize y,z, too: y = (1/w -w)/( C - 4*ln(w) - w^2 + 1/w^2 )^(1/2) z = 2 /( C - 4*ln(w) - w^2 + 1/w^2 )^(1/2) so the curves can be easily sketched. Moreover, we can now compare these to the original system and discover that the relationship between the original 'parameter' x and the current w is dx/dw = -(1/w^2 + 1)/( C - 4*ln(w) - w^2 + 1/w^2 )^(1/2) Actually I've come to be fond of the parameter t = -ln(w); then what we need is dx/dt = (exp(t)+exp(-t))/(C-exp(-2*t)+exp(2*t)+4*t)^(1/2) This doesn't integrate "in finite terms", but whatever the indefinite integral is, we may take for x, and then invert to get t = t(x). This may be composed with the parameterization y = (exp(t)-exp(-t))/(C-exp(-2*t)+exp(2*t)+4*t)^(1/2) z = 2 /(C-exp(-2*t)+exp(2*t)+4*t)^(1/2) to solve the original system. dave ============================================================================== From: jr@redmink.demon.co.uk (John R Ramsden) Subject: Re: Differential Equation Problem Date: Mon, 19 Nov 2001 20:05:43 GMT Newsgroups: sci.math rusin@vesuvius.math.niu.edu (Dave Rusin) wrote: > > The rescalings allow us to assume the problem to solve is just > dy/dx = 1 - y*Sqrt[y^2 + z^2] > dz/dx = - z*Sqrt[y^2 + z^2] In that case taking u, v = y/z, 1/z reduces the system to: du/dx = v dv/dx = sqrt(u^2 + 1) These imply: sqrt(u^2 + 1).du = v.dv so taking: u = (w^2 - 1)/(2.w) one can integrate to obtain: 4.v^2 = w^2 + log (w^4) - 1/w^2 + C Differentiating the preceding (u=) equation WRT x gives: du/dx = (w^2 + 1)/(2.w^2) * dw/dx and equating this to v (expressed in terms of w) yields: / x = | (w^2 + 1) / 2.w^2.sqrt(w^2 + log(w^4) - 1/w^2 + C) dw + D / As u, v are known in terms of w, inverting the latter integral gives u, v and hence y, z in terms of x. Cheers --------------------------------------------------------------------------- John R Ramsden (jr@redmink.demon.co.uk) --------------------------------------------------------------------------- The new is in the old concealed, the old is in the new revealed. St Augustine. ---------------------------------------------------------------------------