MathGroup Archive 2006

[Date Index] [Thread Index] [Author Index]

Search the Archive

Re: Relatively simple, but problematic, non-linear ODE

  • To: mathgroup at smc.vnet.net
  • Subject: [mg72010] Re: Relatively simple, but problematic, non-linear ODE
  • From: ab_def at prontomail.com
  • Date: Fri, 8 Dec 2006 06:17:59 -0500 (EST)
  • References: <el685m$2p9$1@smc.vnet.net>

Alan wrote:
> I am trying to find positive, non-decreasing solutions y(z)
> on the interval (0,z), z > 0 to the ODE:
>
> (y')^2 - rho y y' + (1/4) y^2 = 1,   with y[0] = 0.
>
> where y' = dy/dz,  and the constant parameter |rho| < 1.
> There are generally two solutions. The one I want has the
> same sign as z.
>
> Here is my code:
>
> FSolver[z_, rho_] :=
>   Module[{solnpair, vals, ans},
>       solnpair =  y[x] /. NDSolve[{y'[x]^2 - rho y'[x] y[x] + (1/4) y[x]^2
> == 1,
>           y[0] == 0}, y[x], {x, 0, z}];
>       vals = Re[solnpair /. x -> z];
>       ans = Select[vals, (Sign[#] == Sign[z]) &][[1]];
>    Return[ans];
>
> You can see that y[x] = 2 is also a solution to the ODE.
> The solution I want apparently switches to y = 2
> at some critical value of z. The problem is that this causes Mathematica to
> choke.
> For example, for rho = 0, the exact solution is y(z) = 2 Sin(z/2)
> for |z| < Pi and y(z) = 2 for larger z. I can fix things for rho = 0, since
> I
> am very confident about the answer. But I am less confident about
> the solution for rho not equal to zero.
>
> For example, run my fragment with z = 3,  rho = 1/2 to see the problems.
> Any suggestions on how to get some reliable output
> from NDSolve for this problem will be much appreciated.
>
> Thanks!
> alan

I think if you construct a function that increases from y[0] == 0 to
y[somepoint] == 2 and then stays constant, it will satisfy the
conditions. First find y'[z]:

In[1]:= ode = y'[z]^2 - rho*y[z]*y'[z] + y[z]^2/4 - 1;

In[2]:= Solve[ode == 0, y'[z]]

Out[2]= {{y'[z] -> 1/2*(rho*y[z] - Sqrt[4 - y[z]^2 + rho^2*y[z]^2])},
{y'[z] -> 1/2*(rho*y[z] + Sqrt[4 - y[z]^2 + rho^2*y[z]^2])}}

In[3]:= y'[z] /. % /. y[z] -> 0

Out[3]= {-1, 1}

We want the second solution, because its derivative at z == 0 is
positive. Integrate this new ODE up to the point y[z] == 2:

f[rho_?NonNegative] := Module[{y, z},
  y = y /. First@ NDSolve[
    {y'[z] == 1/2*(rho*y[z] + Sqrt[(rho^2 - 1)*y[z]^2 + 4]),
     y[0] == 0},
    y, {z, 0, Pi},
    StoppingTest -> Re[y[z]] - 2];
  Function @@ {z,
    Piecewise[{{y[z], z < y[Domain[]][[1, 2]]}}, 2]}
]

Verify that this piecewise function satisfies the original ODE:

Plot[ode /. y -> f[rho] /. rho -> .5 // Evaluate,
  {z, 0, 4}, PlotRange -> All]

Unless rho == 0, the solution is continuous but its derivative has a
jump:

Plot[Table[f[rho]'[z], {rho, 0, 1, 1/2}] // Evaluate, {z, 0, 4}]

For -1 < rho < 0 NDSolve already returns a valid solution, which only
approaches the value y[z] == 2 asymptotically.

Maxim Rytin
m.r at inbox.ru


  • Prev by Date: RE: Re: General--Mathematica and Subversion
  • Next by Date: Re: Finding the periphery of a region
  • Previous by thread: Re: Relatively simple, but problematic, non-linear ODE
  • Next by thread: Re: Relatively simple, but problematic, non-linear ODE