MathGroup Archive 2004

[Date Index] [Thread Index] [Author Index]

Search the Archive

NDSolve for Newtonian Orbits Question

  • To: mathgroup at
  • Subject: [mg47297] NDSolve for Newtonian Orbits Question
  • From: "David Park" <djmp at>
  • Date: Fri, 2 Apr 2004 03:31:49 -0500 (EST)
  • Sender: owner-wri-mathgroup at

Dear MathGroup,

I am trying to obtain a numerical solution for Newtonian orbits. (I can solve the
du/dphi equation symbolically with DSolve to obtain the equation of an ellipse, but I also want to know how to do numerical solutions.) I'm having difficulty in knowing how to use NDSolve. For a start I just want to solve for the radius r as a function of time.

Here are the equations. We have an effective potential given by

Veff[M_, h_][r_] = -M/r + (1/2)*(h^2/r^2)

M is the attracting mass and h is the angular momentum. The following is a plot of a particular case.

Plot[Veff[1, 1][r], {r, 0.5, 10},
    PlotRange -> All,
    Frame -> True,
    FrameLabel -> {r, Veff},
    Axes -> False,
    ImageSize -> 450];

The differential equation is given by...

deqn[M_, h_, En_] = En == (1/2)*Derivative[1][r][t]^2 + Veff[M, h][r[t]]
En == h^2/(2*r[t]^2) - M/r[t] + (1/2)*Derivative[1][r][t]^2

where En is the energy. If I pick the energy to be -0.3 and put r[0] somewhere in the orbit range I should obtain a nice elliptical orbit. But when I try to solve I run into all kinds of problems.

NDSolve[{deqn[1, 1, -0.3], r[0] == 1.5}, r, {t, 0, 10}]
r[t_] = r[t] /. Part[%, 2]

Plot[r[t], {t, 0, 5.25}];

It appears that the numerical solution gets stuck at the turning points. How can I obtain an extended solution?

David Park
djmp at 

  • Prev by Date: Re: Abs function
  • Next by Date: Re: How to auto insert Date in footer of each page
  • Previous by thread: Re: Re: Is there any productive way to use Mathematica + pdfLaTeX?
  • Next by thread: Re: NDSolve for Newtonian Orbits Question