Re: Solving Weissinger's ODE
- To: mathgroup at smc.vnet.net
- Subject: [mg104856] Re: Solving Weissinger's ODE
- From: pratip <pratip.chakraborty at gmail.com>
- Date: Thu, 12 Nov 2009 06:01:00 -0500 (EST)
- References: <hde08e$spq$1@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 Nov 11, 10:30 am, Virgil Stokes <v... 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