Solving the path of Earth around Sun

Question

This maybe isn't a universally helpful question. Maybe a little of a code dump. But here goes. I'm trying to solve the path Earth moves around Sun from Earths mass, Suns mass, Earths initial velocity and position. Using NDSolve this set of equations produces a right kind of result:

nSolutionWorks = 
 NDSolve[{
   forceParticle2ExertsOnParticle1[1][t] == (x1[1]'')[t], 
   forceParticle2ExertsOnParticle1[2][t] == (x1[2]'')[t],                    
   forceParticle1ExertsOnParticle2[1][t] == 10 (x2[1]'')[t], 
   forceParticle1ExertsOnParticle2[2][t] == 10 (x2[2]'')[t], 
   forceParticle1ExertsOnParticle2[1][t] == -forceParticle2ExertsOnParticle1[1][t], 
   forceParticle1ExertsOnParticle2[2][t] == -forceParticle2ExertsOnParticle1[2][t], 
   forceParticle2ExertsOnParticle1[1][t] == (10 (-x1[1][t] + x2[1][t]))/(Abs[-x1[1][t] + x2[1][t]]^2 +  Abs[-x1[2][t] + x2[2][t]]^2)^(3/2), 
   forceParticle2ExertsOnParticle1[2][t] == (10 (-x1[2][t] + x2[2][t]))/(Abs[-x1[1][t] + x2[1][t]]^2 +  Abs[-x1[2][t] + x2[2][t]]^2)^(3/2), 
   x1[1][0] == 1, 
   x1[2][0] == 0, 
   x2[1][0] == 0, 
   x2[2][0] == 0, 
   Derivative[1][x1[1]][0] == 0, 
   Derivative[1][x1[2]][0] == 2, 
   Derivative[1][x1[1]][0] + 10 Derivative[1][x2[1]][0] == 0, 
   Derivative[1][x1[2]][0] + 10 Derivative[1][x2[2]][0] == 0}, 
  {x1[1], x1[2], x2[1], x2[2]}, {t, 0, 10}]

The equations use dummy values for the initial values I mentioned. The next equation has these values replaced with real ones. I just can't get it to work though. Any ideas?

nSolutionDoesntWork = NDSolve[
  {forceParticle2ExertsOnParticle1[1][t] == 
    5.9721986`*^24 (x1[1]'')[t], 
   forceParticle2ExertsOnParticle1[2][t] == 
    5.9721986`*^24 (x1[2]'')[t], 
   forceParticle1ExertsOnParticle2[1][t] == 
    1.988435`*^30 (x2[1]'')[t], 
   forceParticle1ExertsOnParticle2[2][t] == 
    1.988435`*^30 (x2[2]'')[t], 
   forceParticle1ExertsOnParticle2[1][
     t] == -forceParticle2ExertsOnParticle1[1][t], 
   forceParticle1ExertsOnParticle2[2][
     t] == -forceParticle2ExertsOnParticle1[2][t], 
   forceParticle2ExertsOnParticle1[1][t] == (
    7.925404384598102`*^44 (-x1[1][t] + 
       x2[1][t]))/(Abs[-x1[1][t] + x2[1][t]]^2 + 
      Abs[-x1[2][t] + x2[2][t]]^2)^(3/2), 
   forceParticle2ExertsOnParticle1[2][t] == (
    7.925404384598102`*^44 (-x1[2][t] + 
       x2[2][t]))/(Abs[-x1[1][t] + x2[1][t]]^2 + 
      Abs[-x1[2][t] + x2[2][t]]^2)^(3/2), 
   x1[1][0] == 1.4956645729619998`*^11, x1[2][0] == 0, x2[1][0] == 0, 
   x2[2][0] == 0, Derivative[1][x1[1]][0] == 0, 
   Derivative[1][x1[2]][0] == 465.10109`, 
   5.9721986`*^24 Derivative[1][x1[1]][0] + 
     1.988435`*^30 Derivative[1][x2[1]][0] == 0, 
   5.9721986`*^24 Derivative[1][x1[2]][0] + 
     1.988435`*^30 Derivative[1][x2[2]][0] == 0},
  {x1[1], x1[2], x2[1], x2[2]},
  {t, 0, 60*60*24*300}
  ]

Evaluating the second NDSolve gives two error messages:

NDSolve::ivres: NDSolve has computed initial values that give a zero residual for the differential-algebraic system, but some components are different from those specified. If you need them to be satisfied, giving initial conditions for all dependent variables and their derivatives is recommended.

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

Answer 1

Eliminate the unnecessary force* variables, which forces NDSolve to use the DAE solver, and just solve it as an ODE. Then you can use higher precision etc., if desired, but in this case you get the same solution using machine precision as higher precision.

newsys =
  ComplexExpand@Eliminate[Rationalize[#, 0] &@{
       forceParticle2ExertsOnParticle1[1][t] == 
        5.9721986`*^24 (x1[1]'')[t], 
       forceParticle2ExertsOnParticle1[2][t] == 
        5.9721986`*^24 (x1[2]'')[t], 
       forceParticle1ExertsOnParticle2[1][t] == 
        1.988435`*^30 (x2[1]'')[t], 
       forceParticle1ExertsOnParticle2[2][t] == 
        1.988435`*^30 (x2[2]'')[t], 
       forceParticle1ExertsOnParticle2[1][
         t] == -forceParticle2ExertsOnParticle1[1][t], 
       forceParticle1ExertsOnParticle2[2][
         t] == -forceParticle2ExertsOnParticle1[2][t], 
       forceParticle2ExertsOnParticle1[1][
         t] == (7.925404384598102`*^44 (-x1[1][t] + 
             x2[1][t]))/(Abs[-x1[1][t] + x2[1][t]]^2 + 
            Abs[-x1[2][t] + x2[2][t]]^2)^(3/2), 
       forceParticle2ExertsOnParticle1[2][
         t] == (7.925404384598102`*^44 (-x1[2][t] + 
             x2[2][t]))/(Abs[-x1[1][t] + x2[1][t]]^2 + 
            Abs[-x1[2][t] + x2[2][t]]^2)^(3/2)},
     {forceParticle2ExertsOnParticle1[1][t], 
      forceParticle2ExertsOnParticle1[2][t], 
      forceParticle1ExertsOnParticle2[1][t], 
      forceParticle1ExertsOnParticle2[2][t]}] /. _Unequal -> True;

nSolutionDoesntWork = NDSolve[{newsys,
   Rationalize[{x1[1][0] == 1.4956645729619998`*^11, x1[2][0] == 0, 
     x2[1][0] == 0, x2[2][0] == 0, Derivative[1][x1[1]][0] == 0, 
     Derivative[1][x1[2]][0] == 465.10109`, 
     5.9721986`*^24 Derivative[1][x1[1]][0] + 
       1.988435`*^30 Derivative[1][x2[1]][0] == 0, 
     5.9721986`*^24 Derivative[1][x1[2]][0] + 
       1.988435`*^30 Derivative[1][x2[2]][0] == 0},
    0]},
  {x1[1], x1[2], x2[1], x2[2]}, {t, 0, 60*60*24*300}
  (*,PrecisionGoal->15, WorkingPrecision->50*)]

ParametricPlot[{{x1[1][t], x1[2][t]}, {x2[1][t], x2[2][t]}} /. 
  First@nSolutionDoesntWork, {t, 0, 60*60*24*300}, AspectRatio -> 1/2,
  PlotPoints -> 1001]

(You can play with the Rationalize commands to see if they are really necessary.)

Answer 2

A somewhat redundant answer to Michael E2's, but providing the simplification.

The problem is not that the numbers are not scaled. (Though it is nice to scale your problems when possible.) The problem is the intermediate forces. You should eliminate them, as well as simplifying the initial conditions. I named your constants for clarity:

NDSolve[{
  x1[1]''[t] == 
   gmm/m1 (x2[1][t] - x1[1][t]) /
    ((x2[1][t] - x1[1][t])^2 + (x2[2][t] - x1[2][t])^2)^(3/2),
  x1[2]''[t] == 
   gmm/m1 (x2[2][t] - x1[2][t]) /
    ((x2[1][t] - x1[1][t])^2 + (x2[2][t] - x1[2][t])^2)^(3/2),
  x2[1]''[t] == -(m1/m2) x1[1]''[t],
  x2[2]''[t] == -(m1/m2) x1[2]''[t],
  x1[1][0] == xi, x1[2][0] == 0,
  x2[1][0] == 0, x2[2][0] == 0,
  x1[1]'[0] == 0, x1[2]'[0] == vi,
  x2[1]'[0] == 0, x2[2]'[0] == -(m1/m2) vi},
 {x1[1], x1[2], x2[1], x2[2]}, {t, 0, tf}]

Then it works fine with both:

m1 = 1;
m2 = 10;
gmm = 10;
xi = 1;
vi = 2;
tf = 2 \[Pi] \[Mu] ((2 \[Mu])/xi - ((1 + m1/m2) vi)^2)^(-3/2) /.
     \[Mu] -> (gmm (m1 + m2))/(m1 m2);

and:

m1 = 5.9721986`*^24;
m2 = 1.988435`*^30;
gmm = 7.925404384598102`*^44;
xi = 1.4956645729619998`*^11;
vi = 465.10109`;
tf = 2 \[Pi] \[Mu] ((2 \[Mu])/xi - ((1 + m1/m2) vi)^2)^(-3/2) /.
     \[Mu] -> (gmm (m1 + m2))/(m1 m2);

Also I set the integration times to the period of the orbit.

For these problems, I would recommend that instead of putting the origin at one of the bodies, that you put the origin at the center of mass.

By the way, your initial velocity for the Earth is extremely small, resulting in a highly elliptical orbit. The Earth almost falls straight down to the Sun. For that case, you should increase the PrecisionGoal for NDSolve[] to get the very near flyby right and to return to close to the same position as where you started. I set it to 13 to get an accurate answer. Machine precision numbers worked fine with that.

I presume that you know that there is a closed form solution for this problem? For two bodies, you don't need to use NDSolve[] at all.