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.