Integrating extremely elliptic orbits displaced from the origin.
- To: mathgroup at smc.vnet.net
- Subject: [mg32536] Integrating extremely elliptic orbits displaced from the origin.
- From: Marco Masi <marco.masi at spiro.fisica.unipd.it>
- Date: Sat, 26 Jan 2002 04:07:53 -0500 (EST)
- Sender: owner-wri-mathgroup at wolfram.com
I tried to integrate numerically normal (2D) Keplerian orbits with
NDSolve. Everything functions fine and well for circular and elliptic
orbits.
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
ellipse).
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
worse.
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
Units*)
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]"}];