MathGroup Archive 2002

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

Search the Archive

Integrating extremely elliptic orbits displaced from the origin.

  • To: mathgroup at
  • Subject: [mg32536] Integrating extremely elliptic orbits displaced from the origin.
  • From: Marco Masi <marco.masi at>
  • Date: Sat, 26 Jan 2002 04:07:53 -0500 (EST)
  • Sender: owner-wri-mathgroup at

I tried to integrate numerically normal (2D) Keplerian orbits with
NDSolve. Everything functions fine and well for circular and elliptic

But I need to integrate extremely elliptic orbits (comets falling towards
the sun with perihelion 0.1 AU and aphelia about 140.000 AU). This
continues to work well as long the orbit isn't displaced from the origin
of coordinates (in the following program for "shift > 100")

It seems that Mathematica is no longer able to integrate orbits correctly
then. The test particle looses or acquires energy (i.e. the orbit spirals
towards the gravity center or away from it instead of becoming an

It seems to me that this isn't because of a lack of WorkingPrecision,
AccuracyGoal, PrecisionGoal, MaxStepsSize or the Method used. I tried many
different parameters and methods but never recovered the ellipse.  
Increasing precision and/or accuracy beyond a threshold makes things even

Does anyone know why this is so? What am I doing wrong?

Thanks to anyone who can help, 
Marco Masi.
The program is as follows:

(*These are various constant*) 
G = 6.67259 10^-11; c = 299792458; M0 = 1.989 10^30;
AU = 149.7 10^9; yr = 3600*24*365 ; Myr = 10^6 * yr; 
pc = c*yr*3.2638; pcau = 206124; 
K = G*(M0*Myr^2)/pc^3; T = 55 Sqrt[K];

(*Shift tells how much the orbit should be displaced from the origin along
the X-axes*)

shift = 100;

(*p is the perihelion and a the aphelion of the orbit, which are
transformed in parsec*)

p = 0.1; a = 140000; perihelion = p*AU/pc; aphelion = a*AU/pc;

(*This solves Newton's classical equations of motion: F = - \nabla
Phi(x,y) /propto - 1 / Sqrt[((x-shift)^2 + y^2)]. Initial velocity x'[0]
is calculated to reach aphelion a from perihelion p*)

orbit = NDSolve[{
x''[t] = -(x[t]-shift)/((x[t]-shift)^2 + y[t]^2)^(3/2),
y''[t] = -y[t]/((x[t]-shift)^2 + y[t]^2)^(3/2), 
x[0] = shift, x'[0] = Sqrt[2*((1/perihelion) - (1/(aphelion + perihelion)),
y[0] = -perihelion, y'[0] = 0}
{x, y}, {t, 0, T}, MaxSteps -> Infinity, AccuracyGoal -> 8,
WorkingPrecision -> 20, PrecisionGoal -> 10];

(*Conversion functions for graphical output from parsec to Astronomical

UAx[t_] := pcau*Part[(x[t]/.orbit),1];
UAy[t_] := pcau*Part[(y[t]/.orbit),1];

ParametricPlot[{UAx[t], UAy[t]}, {t, 0, T}, PlotRange -> All, AspectRatio
-> 1, PlotPoints -> 1000, AxesLabel -> {"X [AU]", "Y [AU]"}];

  • Prev by Date: ReplaceAll doesn't replace
  • Next by Date: Trivial changes to formula
  • Previous by thread: RE: ReplaceAll doesn't replace
  • Next by thread: Trivial changes to formula