Solving Equations of motion; ND Solve stiff system suspected

Question

Im modelling mercury, venus,earth and mars orbiting the sun for a 2 body problem.

Below are the initial position velocity conditions for the planets in question. Distance units are in AU Velocity units are in AU/day

Initial conditions

T = {88.0, 224.7, 365.2, 687.0};(*period in days*)
Ecc = {0.20563069, 0.00677323, 0.01671022, 0.09341233};(*eccentricity of orbit*)   
\[Mu] = 2.9591220823*10^-4;(*standard gravitational parameter of sun*)
M = {174.79439000000002`, 0.`, -2.48284000000001`, 
19.412480000000016`};(*Mean Anomoly*)
\[Psi] = M + 2*Ecc*Sin[M Degree] + 5/4*Ecc^2*Sin[2*M Degree] + 
13/12*Ecc^3*Sin[3*M Degree];(*True Anomoly*)
a = {0.38709893, 0.7233319899999999, 1.00000011, 1.52366231};
r = Table[a[[i]] (1 - Ecc[[i]]^2)/(1 + Ecc [[i]] Cos[\[Psi] [[i]] Degree ]), {i, 4}] (*orbital position*)
rx = Table[r[[i]] Cos[\[Psi][[i]] Degree], {i, 4}](*x component of position*)
ry = Table[r[[i]] Sin[\[Psi][[i]] Degree ], {i, 4}](*y component of position*)
v = Table[Sqrt[\[Mu] (2/r[[i]] - 1/a[[i]])], {i, 4}](*orbital velocity*)
vx = v Cos[\[Psi] Degree](*x component of velocity*)
vy = v Sin[\[Psi] Degree](* y component of velocity*)

Solving using NDSolve

Here i am using NDsolve to solve the ODE's for the equations of motion in x,y direction

planet = Table[NDSolve[{(x''[t] + (\[Mu]*x[t])/(x[t]^2 + y[t]^2)^(
   3/2) == 0, (y''[t] + (\[Mu]*
    y[t])/(x[t]^2 + y[t]^2)^(3/2) == 0, x[0] == rx[[i]], 
 y[0] == ry[[i]], x'[0] == vx[[i]], 
 y'[0] == vy[[i]]}, {x, y}, {t,800, 0.1}], {i,4}];

The problem

If i were to set the intial conditons to the magnitudes of postions and velocity, the simulation works fine, however, its when i use the components of (x,y) for position and velocity that i encounter the error message.

NDSolve::ndsz: At t == 76.8007476284912`, step size is effectively zero; singularity or stiff system suspected.

NDSolve::ndsz: At t == 204.5293849935109`, step size is effectively zero; singularity or stiff system suspected.

NDSolve::ndsz: At t == 333.0376141002715`, step size is effectively zero; singularity or stiff system suspected.

General::stop: Further output of NDSolve::ndsz will be suppressed during this calculation.

What can i do to fix this?

Thank you for reading

Answer 1

It is necessary to correctly calculate the orbital velocity

T = {88.0, 224.7, 365.2, 687.0};(*period in days*)Ecc = {0.20563069, 
  0.00677323, 0.01671022, 
  0.09341233};(*eccentricity of orbit*)\[Mu] = \
0.00029591220828349645438389263285760599`10.122911634152604;\
(*standard gravitational parameter of sun*)M = {174.79439000000002`, 
  0.`, -2.48284000000001`, 
  19.412480000000016`};(*Mean Anomoly*)\[Psi] = 
 M + 2*Ecc*Sin[M Degree] + 5/4*Ecc^2*Sin[2*M Degree] + 
  13/12*Ecc^3*Sin[3*M Degree];(*True Anomoly*)a = {0.38709893, 
  0.7233319899999999, 1.00000011, 1.52366231};
r = Table[
   a[[i]] (1 - Ecc[[i]]^2)/(1 + Ecc[[i]] Cos[\[Psi][[i]] Degree]), {i,
     4}] (*orbital position*);
rx = Table[
   r[[i]] Cos[\[Psi][[i]] Degree], {i, 4}](*x component of position*);
ry = Table[
   r[[i]] Sin[\[Psi][[i]] Degree], {i, 4}](*y component of position*);
v = Table[
   Sqrt[\[Mu] (2/r[[i]] - 1/a[[i]])], {i, 4}](*orbital velocity*);
vx = -v Sin[\[Psi] Degree](*x component of velocity*);
vy = v Cos[\[Psi] Degree](*y component of velocity*);


pl = {Mercury, Venus, Earth, Mars};

planet = Table[
   NDSolve[{x''[t] + (\[Mu]*x[t])/(x[t]^2 + y[t]^2)^(3/2) == 0, 
     y''[t] + (\[Mu]*y[t])/(x[t]^2 + y[t]^2)^(3/2) == 0, 
     x[0] == rx[[i]], y[0] == ry[[i]], x'[0] == vx[[i]], 
     y'[0] == vy[[i]]}, {x, y}, {t, 0, 800}], {i, 4}];

{Table[ParametricPlot[{x[t], y[t]} /. planet[[i]], {t, 0, 800}, 
   PlotLabel -> pl[[i]]], {i, 1, Length[planet]}], 
 ParametricPlot[
  Evaluate[Table[{x[t], y[t]} /. planet[[i]], {i, Length[pl]}]], {t, 
   0, 800}, PlotLegends -> pl]}