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]"}];