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>