How to work with ParametricNDSolve?

Update 3: I got it. A small dumb mistake... It has to be Evaluate[xP[y]'[0] /. s10] like this. Oh well...

Update 2: So if do this Evaluate[xP'[y][0] /. s10] regardless of what you set as y the output is always 0 even though the starting velocity in the direction of x is 1000. Does anybody know what I'm doing wrong? Evaluate[xP[y][t] /. s10] works as it should.

Update: I managed to include variable intial conditions with ParametricNDSolve as user21 suggested. But now the function bVP doesn't work anymore. Now for bVP[y,0] the output is always 1 regardless of y. Also after the swing-by the smaller I set y (the greater the distance), the greater is the output. What needs to be changed to calculate the speed of the probe? I also changed the old code below with the ParametricNDSolve one.

Currently I'm trying to make a simulation of a swing by maneuvre, in which I send a probe around the Moon. For now there is only the probe and the Moon but eventually it will be expanded. I'm doing this using differential equations. The goal of this is to find out how the distance and the starting velocity affects the gain of velocity.

Now I'm having trouble figuring out how to make the initial conditions variable so that I can plot the ratio of inital and final velocity dependent on distance and starting velocity.

This is the working code so far:

grav12x[r1x_, r1y_, m1_, r2x_, r2y_, m2_] :=
G m1 m2 (r2x - r1x)/((r2x - r1x)^2 + (r2y - r1y)^2)^(3/2)
grav12y[r1x_, r1y_, m1_, r2x_, r2y_, m2_] :=
G m1 m2 (r2y - r1y)/((r2x - r1x)^2 + (r2y - r1y)^2)^(3/2)
tmax = 24*60*60*30;
G = 6.67*10^-11;
mMoon = 7.349*10^22;
mProbe = 3000;
s = ParametricNDSolve[{
xM''[t] ==
(grav12x[xM[t], yM[t], mMoon, xP[t], yP[t], mProbe])/mMoon,
yM''[t] ==
(grav12y[xM[t], yM[t], mMoon, xP[t], yP[t], mProbe])/mMoon,
xP''[t] ==
(grav12x[xP[t], yP[t], mProbe, xM[t], yM[t], mMoon])/mProbe,
yP''[t] ==
(grav12y[xP[t], yP[t], mProbe, xM[t], yM[t], mMoon])/mProbe,
xM[0] == 0,
yM[0] == 0,
xP[0] == -100000000,
yP[0] == y,
xM'[0] == -1000,
yM'[0] == 0,
xP'[0] == 1000,
yP'[0] == 0
}, {xM, yM, xP, yP}, {t, 0, tmax},{y}]


Now I want for example yP[0] to be variable so that I can plot this

bVP[y_,t_] :=
Norm[Evaluate[{xP'[y][t], yP[y][t]} /. s]]


which describes the magnitude of the velocity.

I know this is wrong but I'm looking for something like this:

vVerhältnis[k_] := Replace[y, y -> k] && bVP[60*60*24*8]/bVP[0]


I'd like to plot that as a function of, in this case, the vertical distance.

I'll be very thankful if anyone could help me out :)

• Could you please put the code for the constants in your post. For example what are tmax, G. What you are looking for is ParametricNDSolve and then you can make a sensitivity analysis. Nov 30, 2015 at 23:48
• @user21 I added the constants to the code. I'll try it with ParametricNDSolve. Dec 1, 2015 at 15:07
• Welcome to Mathematica.SE! I suggest the following: 0) Browse the common pitfalls question 1) As you receive help, try to give it too, by answering questions in your area of expertise. 2) Read the faq! 3) When you see good questions and answers, vote them up by clicking the gray triangles, because the credibility of the system is based on the reputation gained by users sharing their knowledge. Also, please remember to accept the answer, if any, that solves your problem, by clicking the checkmark sign! Mar 20, 2016 at 17:56

The code is much simpler using vector notation in NDSolve[ ];

tmax   = 5 10^5;
G      = 6.67*10^-11;
mMoon  = 7.349*10^22;
mProbe = 3000;

grav12[r1_, m1_, r2_, m2_] := G m1 m2 (r2 - r1)/((r2 - r1).(r2 - r1))^(3/2)

s = ParametricNDSolve[
{M''[t] == (grav12[M[t], mMoon, P[t], mProbe])/mMoon,
P''[t] == (grav12[P[t], mProbe, M[t], mMoon])/mProbe,
M[0]   == {0, 0},      P[0]  == {-10^9, y},
M'[0]  == {-1000, 0},  P'[0] == {1000, 0}},
{M, P}, {t, 0, tmax}, {y}]

ParametricPlot[(P /. s)[100000][t], {t, 0, 495000}, AspectRatio -> 1]