# Finding the solution for argument of interpolating function for known functional value

Suppose I want to numerically determine the location of an object falling vertically in a gravitational field from an arbitrary height h. I know that I can use NDSolve to get a numerical interpolating function and it's easy enough to plot. However, I want to find the time when the altitude is zero, i.e., when the object hits the surface of the gravitational body.

sol3[hh_] := Module[{h = hh},
NDSolve[{r''[t] == -((gg mm)/(re + r[t])^2), r[0] == h,
r'[0] == 0} /. {re -> 6.37814 10^6, mm -> 5.9742 10^24,
gg -> 6.67430 10^-11}, r, {t, 0, 10000}]
];


I can plot this and see that the object's altitude reaches zero before t = 150 using:

Plot[Evaluate[r[t] /. sol3[10^5]], {t, 0, 150}, PlotRange -> All]


However, I want to know the exact value of t when r[t] equals to zero. I am not sure how to procede as I've tried several approaches but to no avail. Any suggestions would be greatly appreciated.

• "find the time when the altitude is zero" - sounds like a job for WhenEvent[]. Jan 29, 2021 at 17:28
• FindRoot?: mathematica.stackexchange.com/questions/23609/… Jan 29, 2021 at 17:29
• WhenEvent example: mathematica.stackexchange.com/questions/42304/… Jan 29, 2021 at 17:32
• @MichaelE2 - I tried the FindRoot option and kept getting errors like "...is not a list of numbers with dimensions {1} at {x}={0.}..." I haven't tried the WhenEvent route yet but can screw around with that a little. Jan 29, 2021 at 17:37
• For instance, FindRoot[r[t] == 0 /. First@sol3[10^5], {t, 100}] works for me. Jan 29, 2021 at 18:05

@MichaelE2's suggestion was a good one, namely to use WhenEvent to handle the root finding. I couldn't get FindRoot to work properly for a solution that depended upon a user's variable argument (i.e., hh here) but am guessing that was due to my limited capacities with the code.

The following will output the value of t when the object reaches zero altitude:

sol4[hh_] := Module[{h = hh},
Reap[
NDSolve[{r''[t] == -((gg mm)/(re + r[t])^2), r[0] == h,
r'[0] == 0,
WhenEvent[r[t] == 0, Sow[t]]} /. {re -> 6.37814 10^6,
mm -> 5.9742 10^24, gg -> 6.67430 10^-11},
r, {t, 0, 10000}];
][[2, 1]]
];


So when I enter 100 km for hh I get a time of 144.711 seconds, which seems reasonable since we ignore air resistance.

• Additionally, ParametricNDSolve[] is convenient for a situation like yours: With[{re = 6.37814*^6, mm = 5.9742*^24, gg = 6.67430*^-11}, sol3 = ParametricNDSolveValue[{r''[t] == -((gg mm)/(re + r[t])^2), r[0] == h, r'[0] == 0, WhenEvent[r[t] == 0, Sow[t]; "StopIntegration", "DetectionMethod" -> "Interpolation", "LocationMethod" -> "Brent"]}, r, {t, 0, 10000}, h]]. Then try Reap[sol3[100000]] Jan 29, 2021 at 17:59
• pos[t_] = r[t] /. sol3[10^5]; FindRoot[pos[t] == 0, {t, 140}] This gives: {t -> 144.711} Jan 29, 2021 at 19:28

" I couldn't get FindRoot to work properly ..." Try this

rsol[h_?NumericQ] :=
r /. First@
NDSolve[{r''[t] == -((gg mm)/(re + r[t])^2), r[0] == h,
r'[0] == 0} /. {re -> 6.37814 10^6, mm -> 5.9742 10^24,
gg -> 6.67430 10^-11}, r, {t, 0, 10000}] // Quiet

Plot[Evaluate[rsol[10^5][t]], {t, 0, 150}, PlotRange -> All]

FindRoot[rsol[10^5][t] == 0, {t, 140}]

(*   {t -> 144.711}   *)


Edit you can also look for h at given t

FindRoot[rsol[h][144.7108864471939] == 0, {h, 90000},
AccuracyGoal -> 6]

(*   {h -> 100000.}   *)
`
• Thanks! I really wish I had more time to play with Mathematica. I often have several month droughts while I'm doing my job and then get a few moments every so often to come back and play with it. So my memory fails me and I feel like I have to re-learn things. Jan 29, 2021 at 18:23