Finding Minimium Periapse Distance Between Spacecraft and Planet

Question

I have the following simple model that simulates a spacecraft orbiting Mars. What I'm hoping to do is be able to find the periapse radius of the spacecraft as I did below, except instead of using tables I'll use the interpolating functions that are given as output from NDSolve (would using interpolating functions to find the periapse radius be faster than using tables?):

G = 6.672*10^-11;
m[1] = 6.4185*10^23;
m[2] = 100;
p[1] = {0, 0};
p[2] = {1000000, 1000000};
v[1] = {0, 2500};
v[2] = {0, 0};
tmax = 1000;


soln = NDSolve[{
    x[1]''[t] == -(G m[1] (x[1][t] - x[2][t]))/((x[1][t] - x[2][t])^2 + (y[1][t] - y[2][t])^2)^(3/2),
    y[1]''[t] == -(G m[1] (y[1][t] - y[2][t]))/((x[1][t] - x[2][t])^2 + (y[1][t] - y[2][t])^2)^(3/2),
    x[2]''[t] == -(G m[2] (x[2][t] - x[1][t]))/((x[2][t] - x[1][t])^2 + (y[2][t] - y[1][t])^2)^(3/2),
    y[2]''[t] == -(G m[2] (y[2][t] - y[1][t]))/((x[2][t] - x[1][t])^2 + (y[2][t] - y[1][t])^2)^(3/2),
    x[1][0] == p[1][[1]], y[1][0] == p[1][[2]], x[2][0] == p[2][[1]], 
    y[2][0] == p[2][[2]], x[1]'[0] == v[1][[1]], 
    y[1]'[0] == v[1][[2]], x[2]'[0] == v[2][[1]], 
    y[2]'[0] == v[2][[2]]}, {x[1][t], y[1][t], x[2][t], y[2][t], 
    x[3][t], y[3][t], x[4][t], y[4][t]}, {t, 0, tmax}, 
   Method -> "StiffnessSwitching", AccuracyGoal -> 18, 
   PrecisionGoal -> 18, MaxSteps -> 10000000];

ParametricPlot[{{{x[1][t], y[1][t], {x[2][t], y[2][t]}}}} /. soln, {t,
    0, tmax}, Prolog -> {Red, Disk[{p[2][[1]], p[2][[2]]}, 50000]}, 
  AxesLabel -> {x, y}, ImageSize -> Large]

Animate[ParametricPlot[{{x[1][t], y[1][t]}, {x[2][t], y[2][t]}} /. 
     soln /. t -> a, {t, Max[0, a - 1000], a}, 
   Prolog -> {Red, Disk[{p[2][[1]], p[2][[2]]}, 50000]}, 
   AxesLabel -> {x, y}, Axes -> False, ImageSize -> Large], {a, 0, 
   tmax}, AnimationRate -> 10]


   dt = 1;
MarsPosition = Table[{x[1][t], y[1][t]} /. soln, {t, 0, tmax, dt}] ;
SpaceCraftPosition = Table[{x[2][t], y[2][t]} /. soln, {t, 0, tmax, dt}];
dxy = Sqrt[(MarsPosition - SpaceCraftPosition)^2];
dr = Table[Norm[dxy[[i]]], {i, 1, Length[dxy]}];
mindr = Min[dr] (*Find closest approach of spacecraft to Mars*)
mindrindex = Flatten[Position[dr, mindr]][[1]] (*Finds index position of mindr in dr*)

I've already had a look at a few topics that used FindMinValue, FindMinimum and NMinimize (Finding the Minimum value of an interpolating function, Find minimum value of function which calls NDsolve, Find maximum value of interpolation function - obviously wrong result), but could not get them to work as I would always get syntax errors that I could not solve. Would using these functions be the correct way of achieving what I'd like to do, and if so, could anyone help me figure out how to get past the syntax errors? Any help would be brilliant.

Answer 1

G = 6672*10^-14;
m[1] = 64185*10^19;
m[2] = 100;
p[1] = {0, 0};
p[2] = {10^6, 10^6};
v[1] = {0, 2500};
v[2] = {0, 0};
tmax = 1000;

soln = NDSolve[{
    x[1]''[t] == -(G m[1] (x[1][t] - x[2][t]))/
      Norm[{x[1][t], y[1][t]} - {x[2][t], y[2][t]}]^3,
    y[1]''[t] == -(G m[1] (y[1][t] - y[2][t]))/
      Norm[{x[1][t], y[1][t]} - {x[2][t], y[2][t]}]^3,
    x[2]''[t] == -(G m[2] (x[2][t] - x[1][t]))/
      Norm[{x[1][t], y[1][t]} - {x[2][t], y[2][t]}]^3,
    y[2]''[t] == -(G m[2] (y[2][t] - y[1][t]))/
      Norm[{x[1][t], y[1][t]} - {x[2][t], y[2][t]}]^3,
    x[1][0] == p[1][[1]], y[1][0] == p[1][[2]],
    x[2][0] == p[2][[1]], y[2][0] == p[2][[2]],
    x[1]'[0] == v[1][[1]], y[1]'[0] == v[1][[2]],
    x[2]'[0] == v[2][[1]], y[2]'[0] == v[2][[2]]},
   {x[1][t], y[1][t], x[2][t], y[2][t],
    x[3][t], y[3][t], x[4][t], y[4][t]},
   {t, 0, tmax},
   Method -> "StiffnessSwitching",
   AccuracyGoal -> 18,
   PrecisionGoal -> 18,
   MaxSteps -> 10000000];

mindr = NMinimize[Join[
   Norm[{x[1][t], y[1][t]} - {x[2][t], y[2][t]}] /. soln,
   {0 <= t <= tmax}], t]

{76703.6, {t -> 909.167}}

ParametricPlot[
 {{x[1][t], y[1][t]}, {x[2][t], y[2][t]}} /. soln,
 {t, 0, tmax},
 Prolog -> {
   Red, Disk[{p[2][[1]], p[2][[2]]}, 50000],
   Magenta, AbsolutePointSize[6],
   Point[{x[1][t], y[1][t]} /. soln /. mindr[[-1]]]},
 AxesLabel -> {x, y},
 ImageSize -> Large]

Answer 2

What you want to do is called "event location" and is realized with NDSolve using WhenEvent. In principle you give it a predicate that is true when the spaceship is at periapsis and NDSolve uses a root finding method to figure out exactly when this happens.

G = 6.672*10^-11;
m[1] = 6.4185*10^23;
m[2] = 100;
p[1] = {0, 0};
p[2] = {1000000, 1000000};
v[1] = {0, 2500};
v[2] = {0, 0};
tmax = 650;

system = {
   {x[1]''[t], y[1]''[t]} == -G m[1] {x[1][t] - x[2][t], y[1][t] - y[2][t]}/Norm[{x[1][t] - x[2][t], y[1][t] - y[2][t]}]^3,
   {x[2]''[t], y[2]''[t]} == -G m[2] {x[2][t] - x[1][t], y[2][t] - y[1][t]}/Norm[{x[1][t] - x[2][t], y[1][t] - y[2][t]}]^3,
   {x[1][0], y[1][0]} == p[1],
   {x[2][0], y[2][0]} == p[2],
   {x[1]'[0], y[1]'[0]} == v[1],
   {x[2]'[0], y[2]'[0]} == v[2],
   WhenEvent[
    {x[1][t] - x[2][t], y[1][t] - y[2][t]}.{x[1]'[t] - x[2]'[t], y[1]'[t] - y[2]'[t]} == 0,
    Sow[t]
    ]
   };

initial = {
   x[1][t],
   y[1][t],
   x[2][t],
   y[2][t],
   x[3][t],
   y[3][t],
   x[4][t],
   y[4][t]
   };

{soln, periapsis} = Reap@NDSolve[
    system,
    initial,
    {t, 0, tmax},
    Method -> "StiffnessSwitching",
    AccuracyGoal -> 18,
    PrecisionGoal -> 18,
    MaxSteps -> 10000000
    ];

pts = Table[Flatten[{x[1][t], y[1][t]} /. soln], {t, Flatten@periapsis}]

ParametricPlot[
 {{x[1][t], y[1][t], {x[2][t], y[2][t]}}} /. soln,
 {t, 0, tmax},
 Prolog -> {
   Red,
   Disk[p[2], 50000],
   Green,
   PointSize[Large],
   Point[pts]
   },
 AxesLabel -> {x, y}
 ]

The predicate that I'm using is when the dot product between the velocity vector and the position vector changes sign. This produces two candidates; the periapsis and the apoapsis, as can be seen in the figure.

Note: I don't know which root finding algorithm NDSolve uses but conceptually it will be similar to the secant method. Therefore if you request a high precision is might tell you that it couldn't find the exact place at the precision that you requested because it required too many iterations. But it will tell you how well it can approximate it.