2
$\begingroup$

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.

Martian periapse

$\endgroup$
  • 2
    $\begingroup$ I think your system would benefit greatly from a judicious rescaling of your variables. $\endgroup$ – J. M. is away Jun 13 '15 at 13:13
4
$\begingroup$
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]

enter image description here

$\endgroup$
  • $\begingroup$ Brilliant, thank you Bob, worked like a charm. I'd also like to be able to get the x and y velocity components at periapse, but how would I isolate 909.167 from {t -> 909.167}? I tried mindr[[2]] which produced {t -> 909.167}, so x[1]'[t] /. soln /. t ->mindr[[2]] unfortunately did not work. $\endgroup$ – indigoblue Jun 13 '15 at 18:12
  • $\begingroup$ Nevermind, using Last[Last[mindr[[2]]]] did the trick. $\endgroup$ – indigoblue Jun 13 '15 at 20:34
4
$\begingroup$

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

Mathematica graphics

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.

$\endgroup$
  • 1
    $\begingroup$ The built-in root-finder is in fact Brent's method, which is effectively a hybrid of bisection, the secant method, and inverse quadratic interpolation. $\endgroup$ – J. M. is away Jun 13 '15 at 15:01

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.