0
$\begingroup$

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.

$\endgroup$
7
  • 1
    $\begingroup$ "find the time when the altitude is zero" - sounds like a job for WhenEvent[]. $\endgroup$ Jan 29, 2021 at 17:28
  • $\begingroup$ FindRoot?: mathematica.stackexchange.com/questions/23609/… $\endgroup$
    – Michael E2
    Jan 29, 2021 at 17:29
  • 1
    $\begingroup$ WhenEvent example: mathematica.stackexchange.com/questions/42304/… $\endgroup$
    – Michael E2
    Jan 29, 2021 at 17:32
  • $\begingroup$ @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. $\endgroup$ Jan 29, 2021 at 17:37
  • 1
    $\begingroup$ For instance, FindRoot[r[t] == 0 /. First@sol3[10^5], {t, 100}] works for me. $\endgroup$
    – Michael E2
    Jan 29, 2021 at 18:05

2 Answers 2

1
$\begingroup$

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

$\endgroup$
2
  • 1
    $\begingroup$ 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]] $\endgroup$ Jan 29, 2021 at 17:59
  • 1
    $\begingroup$ pos[t_] = r[t] /. sol3[10^5]; FindRoot[pos[t] == 0, {t, 140}] This gives: {t -> 144.711} $\endgroup$ Jan 29, 2021 at 19:28
1
$\begingroup$

" 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.}   *)
$\endgroup$
1
  • $\begingroup$ 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. $\endgroup$ Jan 29, 2021 at 18:23

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.