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

Question

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.

I start with something like the following for Earth:

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[]. — J. M.'s missing motivation, Jan 29 '21 at 17:28
FindRoot?: https://mathematica.stackexchange.com/questions/23609/using-findroot-with-an-interpolating-function — Michael E2, Jan 29 '21 at 17:29
WhenEvent example: https://mathematica.stackexchange.com/questions/42304/finding-the-intersection-of-a-curve-with-an-interpolation-function/42310#42310 — Michael E2, Jan 29 '21 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. — honeste_vivere, Jan 29 '21 at 17:37
You'd have to show what you tried, if you want someone to fix the code. You might be missing a First somewhere. — Michael E2, Jan 29 '21 at 18:01
For instance, FindRoot[r[t] == 0 /. First@sol3[10^5], {t, 100}] works for me. — Michael E2, Jan 29 '21 at 18:05

score 1 · Answer 1 · answered Jan 29 '21 at 17:45

1

@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.

answered Jan 29 '21 at 17:45

honeste_vivere

243
3
14

1

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]] – J. M.'s missing motivation Jan 29 '21 at 17:59
1

pos[t_] = r[t] /. sol3[10^5]; FindRoot[pos[t] == 0, {t, 140}] This gives: {t -> 144.711} – Daniel Huber Jan 29 '21 at 19:28

Akku14 · Answer 2 · 2021-01-29T18:55:07.967

1

" 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.}   *)

edited Jan 29 '21 at 18:55

answered Jan 29 '21 at 18:17

Akku14

17,287
14
32

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. – honeste_vivere Jan 29 '21 at 18:23

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

2 Answers2