Approximate with an absolute error less than $0.01$s, the necessary time for the object to reach the ground

Question

An object that falls freely through the air is subject not only to the gravitational force but also to (force of) air resistance. If the object has mass $m$ and is let to fall from height $h_0$, after $t$ seconds, the height of the object is:
$$h(t)=h_0-\frac{mg}{k}t+\frac{m^2g}{k^2}(1-e^\frac{-kt}{m})$$

where $g = 9.8\frac{m}{s^2}$ is the gravitational acceleration and $k$ is the air resistance coefficient.
$h_0 = 300$ m
$m = 0.2$ kg
$k = 0.05 \frac{kg}{s}$
Approximate with an absolute error less than $0.01$s, the necessary time for the object to reach the ground. For the solution, use a suitable numerical algorithm to implement in Maple. Then solve the problem using the predefined function corresponding from Maple and compare the results.

To be fair, I don't even know how to begin. This problem looks like physics but I got it from my Numerical Analysis course; the problem is translated so sorry if there are some poorly worded things in there, tried my best.

Thanks in advance

Answer 1

This is an exercise in numeric root-finding.

I am going to assume that you are new to Maple, and also that your numerical analysis requirements here is on the level of an introductory course, involving the most basic algorithms and whatever calculus of other math lies behind them.

You can set up the function, then plot it to see roughly where it attains zero height (crosses the x-axis), call the stock Maple numeric solver (fsolve) on it, then write you own simple iterative root-finder.

Below, I choose the bisection method, and used the plot to see where I could find initial end-points at which the function had opposite signs (ie. negative and positive).

Here's the function. I don't really need to put all those units in, and I have them all resolve away. (It's a bit of a sanity check on my work. But you will check all my work anyway, right?)

restart;
form := h0 - m*g/k*t +m^2*g/k^2*(1-exp(-k*t/m)):
this := eval(form,[h0 = 300*Unit(m),
                   m = 0.2*Unit(kg),
                   k = 0.05*Unit(kg/s),
                   g = 9.8*Unit(m/s^2),
                   t = t*Unit(s)]):
h := unapply(simplify(this/Unit(m)), t);

      h := t -> 456.8000000-39.20000000*t
                -156.8000000*exp(-.25*t)

Now plot it, to get the lay of the land,

plot(h(t), t=0..15);

Now use the Maple command fsolve, to find an approximate root.

fsol := fsolve(h(t), t=0..12);

          fsol := 11.42301119

For the bisection method we need two initial endpoints at which this continuous h has opposite signs, which ensures that there is a root between them. Plotting can help with that.

h(0), h(15);

       300.0000000, -134.8875826

So we can use a=0 and b=15 in a bisection method, since h(0) is of the opposite sign than h(15).

Now you just have to code a simple bisection procedure. As soon an the "latest" endpoints are less than 0.01 apart then you are guaranteed that the midpoint is within 0.01 of an actual root.

This is a simple-minded implementation, which you should try hard to understand. The only bells&whistles are the printf to show how it's proceeding. You can ignore that. But you really should try and understand the flow of the procedure, and that's part of what the exercise is supposed to teach you. Also, check for mistakes!

bis := proc(H, aval::numeric, bval::numeric)
  local a, b, i, r, w;
  a := aval; b := bval;
  for i from 1 to 20 do
    r := evalf((a+b)/2);
    w := evalf(abs(b-a));
    printf("iteration %ld: midpoint = %.4f  width = %.4f\n",
           i,r,w);
    if w < 0.01 then return r; end if;
    if signum(H(r)) = signum(H(a)) then
      a := r;
    else
      b := r;
    end if;
  end do;
  return FAIL;
end proc:

Now we can run it, passing the bis procedure our operator h and the two initial end-points.

bsol := bis(h, 0, 15);

iteration 1: midpoint = 7.5000  width = 15.0000
iteration 2: midpoint = 11.2500  width = 7.5000
iteration 3: midpoint = 13.1250  width = 3.7500
iteration 4: midpoint = 12.1875  width = 1.8750
iteration 5: midpoint = 11.7188  width = 0.9375
iteration 6: midpoint = 11.4844  width = 0.4688
iteration 7: midpoint = 11.3672  width = 0.2344
iteration 8: midpoint = 11.4258  width = 0.1172
iteration 9: midpoint = 11.3965  width = 0.0586
iteration 10: midpoint = 11.4111  width = 0.0293
iteration 11: midpoint = 11.4185  width = 0.0146
iteration 12: midpoint = 11.4221  width = 0.0073

         bsol := 11.42211914

So the value now assigned to bsol should be within 0.01 of a root. Compare with how much it differs from the fsolve result. Ask yourself how many digits to the right of the decimal point are correct for bsol, etc.

Answer 2

Hint.

When $h(t^*) = 0$

$$ e^{-\frac{k t^*}{m}} = 1-\left(\frac{m g t^*}{k}-h_0\right)\frac{k^2}{m^2 g} = 1-\frac{k t^*}{m}+\frac{1}{2!}\left(\frac{k t^*}{m}\right)^2+\cdots.. $$

Note that $e^{-x} = 1 - x + \frac{1}{2!}x^2-\frac{1}{3!}x^3+\cdots +$ develops in a series which converges absolutely, so the approximation obtained by truncation can be easily estimated...

Answer 3

You know that $h(0)=h_0>0$. If one cancels the first three terms, leading to $$ t_1=\frac{kh_0}{mg}+\frac{m}{k}, $$ then only the last negative term remains, $h(t_1)=-\frac{m^2g}{k^2}e^{-kt_1/m}$.

This means that $[0,t_1]$ is a root-bracketing interval, you can now reduce this subsequently using bisection, regula falsi or any bracketing method.

Then, as example, the Illinois variant of the regula falsi method produces an iteration sequence that starts as

       point        function value
a:  0.000000000 ->         300
b: 11.653061224 ->    -8.51395   
c: 11.331475889 ->     3.37944  
c: 11.422852547 ->  0.00586098