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