8
$\begingroup$

I have the following code that is used to simulate a spacecraft orbiting Mars, and uses ParametricNDSolve to vary the spacecraft's initial conditions using a parameter value p:

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

soln = ParametricNDSolve[{
    xm''[t] == -(G m[2] (xm[t] - xsc[t]))/((xm[t] - xsc[t])^2 + (ym[t] - ysc[t])^2)^(3/2),
    ym''[t] == -(G m[2] (ym[t] - ysc[t]))/((xm[t] - xsc[t])^2 + (ym[t] - ysc[t])^2)^(3/2),
    xsc''[t] == -(G m[1] (xsc[t] - xm[t]))/((xsc[t] - xm[t])^2 + (ysc[t] - ym[t])^2)^(3/2),
    ysc''[t] == -(G m[1] (ysc[t] - ym[t]))/((xsc[t] - xm[t])^2 + (ysc[t] - ym[t])^2)^(3/2),

    xm[0] == p[1][[1]], ym[0] == p[1][[2]], xsc[0] == p[2][[1]], 
    ysc[0] == p[2][[2]], xm'[0] == v[1][[1]], ym'[0] == v[1][[2]], 
    xsc'[0] == v[2][[1]] + 1000 p, ysc'[0] == v[2][[2]] + 1000 p},

   {xm[t], ym[t], xsc[t], ysc[t]}, {t, 0, tmax}, {p}, 
   Method -> "StiffnessSwitching", AccuracyGoal -> 18, 
   PrecisionGoal -> 18, MaxSteps -> 10000000];

   ParametricEscapeTrajectories[t_] = ParametricPlot[Evaluate[Table[{xsc[t][p], ysc[t][p]} /. soln, {p, 0, 5,1}]], {t, 0, tmax}, Prolog -> {Red, Disk[{p[1][[1]], p[1][[2]]}, 50000]}]

Below is an ensemble of results created by varying the parameter p from 0 to 5 as shown by the ParametricEscapeTrajectories function.

enter image description here

What I'm hoping is to be able to do something similar to what was done at this link, where the NMinimize function was used to find the closest approach distance between Mars and the spacecraft (the so-called periapse radius). However, this time I'm hoping to be able to find the value of closest approach, time of closest approach AND the value of the parameter p that gives me the closest approach out of all possible parameter values between 0 and 5. So, for example, it may look something like this:

{71729.9, {t -> 118.096}, {p -> 5}}

So far I've tried things such as the following, but unfortunately have not been successful:

NMinimize[Join[Norm[{xm[t][p], ym[t][p]} - {xsc[t][p], ysc[t][p]}] /. soln, {t, 0, tmax}, {p, 0, 5, 1}], {t, p}]

Any help would be greatly appreciated, thanks very much.

$\endgroup$
9
$\begingroup$
NMinimize[{Norm[{xm[t][p], ym[t][p]} - {xsc[t][p], ysc[t][p]}] /.  soln, 
           0 < t < tmax && 0 < p < 5}, 
           {t, p}, Method -> "SimulatedAnnealing"]

(* {71729.9, {t -> 118.095, p -> 5.}} *)

Edit

Please note that your system is more elegant and easier to understand when written in vector form:

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

soln = ParametricNDSolveValue[{
   posM''[t]  == -G m[2] (posM[t] - posSc[t])/  Tr[(posM[t] - posSc[t])^2]^(3/2),
   posSc''[t] == - m[1]/m[2] posM''[t],
   posM[0]  == p[1],  posSc[0] == p[2],
   posM'[0] == v[1],  posSc'[0] == v[2] + 1000 r },
   {posM, posSc}, {t, 0, tmax}, {r}, 
   Method -> "StiffnessSwitching", AccuracyGoal -> 18, 
   PrecisionGoal -> 18, MaxSteps -> 10000000]

ParametricPlot[Evaluate[Table[Last@Through[soln[p][t]], {p, 0, 5, 1}]], {t, 0, tmax}]

Mathematica graphics

$\endgroup$
  • $\begingroup$ Brilliant, thanks very much belisarius. May I ask why simulated annealing is used instead of letting NMinimize choose its own method? $\endgroup$ – indigoblue Jul 11 '15 at 17:03
  • $\begingroup$ @indigoblue I found it faster ... $\endgroup$ – Dr. belisarius Jul 11 '15 at 17:09

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.