Re: Solving Weissinger's ODE
- To: mathgroup at smc.vnet.net
- Subject: [mg104857] Re: [mg104842] Solving Weissinger's ODE
- From: Pratip Chakraborty <pratip.chakraborty at gmail.com>
- Date: Thu, 12 Nov 2009 06:01:11 -0500 (EST)
- References: <200911110929.EAA29439@smc.vnet.net>
Hi Virgil, As far as I understand with respect to your initial value problem
NDSolve is behaving completely in a correct manner. First thing you may not
have noticed is the fact that the left hand side of the ODE can be
considered as a cubic polynomial in term of y'[t]. As expected Mathematica
first factorize the LHS of the ODE and comes up with three linear factors. t
* (y[t])^2 * (y'[t])^3 - (y[t])^3 * (y'[t])^2 + t * (t^2 + 1) * y'[t] - t^2
*y[t] ==0 =>(y'[t]-a1)(y'[t]-a2)(y'[t]-a3)==0 Now as one can guess one root
among the {a1,a2,a3} will be real and the other two will be complex
conjugates. You can see this very fact if we try to find the y'[t] at t==1.
Actually this y'[1] will be used in the first iteration step by the finite
difference scheme (for example explicit Euler) of NDSolve to find
y[1+h]=y[1]+y'[1]*h on the solution path y[t]. But if we look at the y'[1]
we can find that there exists three distinct values of it at t=1 (one real
and other two complex). This simply implies that there will be three
branches of numerical solutions for y[t] corresponding to different values
of y'[1]. You may evaluate the following to see the above fact:
Solve[t*(y[t])^2*(y'[t])^3 - (y[t])^3*(y'[t])^2 + t*(t^2 + 1)*y'[t] -
t^2*y[t] == 0, y'[t]] /. t -> 1 /. y[1] -> Sqrt[3/2] // N At last I must say
that as you already knew the analytic-real solution there could be another
way to check if any of the numerical solution produced by NDSolve were
meaningful or not. If you plot the three branches of the solution you will
see only one of them can be plotted in the real Cartesian plane. The other
two complex solution will produce empty graphs. You may evaluate the
following to see the above fact:
s=NDSolve[{t*(y[t])^2*(y'[t])^3-(y[t])^3*(y'[t])^2+t*(t^2+1)*y'[t]-t^2*y[t]==0,y[1]==Sqrt[3/2]},y,{t,1,10}];
pic=Table[Plot[y[t]/.s[[i]],{t,1,10},PlotRange-> All,Frame-> True],{i,1,3}]
In our case only the first solution is the real solution you are looking
for. This is simply proved once we plot the first numerical solution along
with the analytic-real solution you provided. You may evaluate the following
to see the above fact: picExact=Plot[Sqrt[t^2+1],{t,1,10},PlotRange->
All,Frame-> True,PlotStyle-> {Directive[Red,Dashed,Thick]}];
Show[pic[[1]],picExact] Last but not least is the fact that one can verify
is whether the other two numerical solutions of the given dynamical system
are really complex conjugates or not. We can witness this just by plotting
the two solutions on the complex plane. One can see that they occur as
"Mirror Reflection" of each other with respect to the real axis. You may
evaluate the following to see the above fact:
Show[ParametricPlot[Evaluate[{Re[y[t]],Im[y[t]]}/.s[[2]]],{t,1,10},PlotStyle->{{Thick,Red}},PlotRange->All],ParametricPlot[Evaluate[{Re[y[t]],Im[y[t]]}/.s[[3]]],{t,1,10},PlotStyle->{{Thick,Blue}},PlotRange->All],Frame->
True,PlotRange->All] I hope this answers your question. Regards, Pratip
On Wed, Nov 11, 2009 at 10:29 AM, Virgil Stokes <vs at it.uu.se> wrote:
> I can not see why the following does not work as expected,
>
> s = NDSolve[{t * (y[t])^2 * (y'[t])^3 - (y[t])^3 * (y'[t])^2 + t *
> (t^2 + 1) * y'[t] - t^2 *y[t] == 0, y[1] == Sqrt[3/2]}, y[t], {t, 1, 10}]
>
> Note, the solution to this nonlinear, non-autonomous, implicit ODE for
> initial condition y[1] = Sqrt[3/2] is just y[t] = Sqrt[t^2 + 1].
>
> Any suggestions on how to obtain the solution (either analytic or
> numerical) would be appreciated.
>
> --V. Stokes
>
>
>
- References:
- Solving Weissinger's ODE
- From: Virgil Stokes <vs@it.uu.se>
- Solving Weissinger's ODE