ri = 5; h = 3.33; mu = 5.678; tmax = 5.;
eqn = {R''[t] - R[t] TH'[t]^2 == -mu/R[t]^2, R'[0] == 4.5, R[0] == ri,
    TH'[t] R[t]^2 == h, TH[0] == 3 Pi/4};
NDSolve[eqn, {TH, R}, {t, 0, tmax}];
{r[u_], th[u_]} = {R[u], TH[u]} /. First[%];
ParametricPlot[r[t] {Cos[th[t]], Sin[th[t]] }, {t, 0, tmax}]
Plot[{r[t], th[t]}, {t, 0, tmax}]

Tried to model planet motion using coupled $ (r,\theta) $ standard Newton equations with some constants.

$$ \ddot {r} - \dot{\theta}^2\,r = - \mu/r^2 ;\, r^2 \, \dot{\theta} = h \,;$$

Results do not look like any conic. $\theta $ does not change, where may be the error? Thanks for finding the error. Should boundary condition for $ \dot r$ be from its integrand?

EDIT 1 :

It works alright with a second order $\theta$. There is no coding error.

$$ \ddot {r} - \dot{\theta}^2\,r = - \mu/r^2 ;\quad 2 \dot {r} {\dot \theta} + \ddot {\theta} r =0;\, $$

eqn={R''[t]-R[t] TH'[t]^2==-mu/R[t]^2, R[t] TH''[t]+2R'[t] TH'[t]==0,R'[0]==-.5,R[0]==ri,TH'[0]==h/ri^2,TH[0]==3Pi/4}
Plot[th'[t] r[t]^2,{t,0,tmax}]

I am at a total loss to understand why Runge-Kutta integration needs $\theta$ to match $r$ to a second order for it to work. Shall appreciate any comments.


As mentioned in comments there was no error at all , neither in code nor in the order of $\theta$ in ODE. Initial value R'[0]= 4.5 gives a hyperbola in its asymptote region that looks like a straight radial ray. Next when initial condition as R'[0]= - 0.5 was given resulting in an ellipse, the earlier non-recognition (by all!) could be recognized.The question is resolved now and can be closed.


So it was due to shape recognition error as can be seen by varying initial $r^{\prime}$ and sign of $t_{max}$. The following are integrands of single order $\theta$ ODEs. Numerically the first orbit has $ e= 0.6094 $

Orbits with $ r^{\prime}s $  varied]1

  • $\begingroup$ @Narasimham Your syntax is fine, I made a mistake. Sorry for the confusion. $\endgroup$ – Jack LaVigne Aug 25 '16 at 13:24
  • $\begingroup$ OK.. I am wondering about the need for second order derivative of $\theta$, I remember it worked before with first order. $\endgroup$ – Narasimham Aug 25 '16 at 13:27
  • $\begingroup$ Got it now!! $R'[0] = - 0.5 $ gives the same solution even for ODE as the first order $ \theta $. Also $R'[0] = 4.5 $ is straight portion of hyperbola asymptote that I mistook for a straight line. $\endgroup$ – Narasimham Aug 25 '16 at 14:12
  • $\begingroup$ For me both versions (first or second order in $\theta$) run and one gets identical results with identical inital conditions. What was your problem with the first oder $\theta$-System? $\endgroup$ – N0va Aug 25 '16 at 14:32
  • $\begingroup$ I have no problem remaining now. As explained in comments above I thought the somewhat straight line orbit was wrong as it could not be recognized as part of hyperbola orbit . Now the changed boundary conditions result in a recognizable ellipse. $\endgroup$ – Narasimham Aug 25 '16 at 15:25