Differential Equations in the New Millennium:The Parachute Problem

Douglas B. Meade (meade@math.sc.edu)
Department of Mathematics
University of South Carolina
Columbia, SC 29229
Allan A. Struthers (struther@math.mtu.edu)
Department of Mathematical Sciences
Michigan Technical University
Houghton, MI 49931

This notebook is a CAS supplement to a paper submitted to a special Mathematics
Education of Engineers issue of the International Journal of Engineering Education
(co-edited by Prof. Shirley Pomeranz, Univ.of Tulsa).  A parallel Maple worksheet
is also available.

Introduction

The position (height above ground level), velocity, acceleration, and jerk (rate of change of acceleration) provide a good description of the motion of the skydiver.  Numerical solutions of the various initial value problems considered are computed using the built in Mathematica numerical differential equation solver NDSolve.  This solver is a general purpose variable stepsize (by default the stepsize is adjusted to try to give 16 place accuracy)  numerical solver based on the LSODE package. NDSolve returns InterpolatingFunction objects: with the default settings these are simple cubic splines through the sampled points.  Note, the Mathematica plot routine recursively resamples (the deafault intial sample size is 25 with a default maximum depth of 30 recursive resamplings) where functions have high curvature: as a result it typically resolves a simple discontinuity well representing it by what appears to be a vertical line.

Linear and Quadratic Models:
Analytical Solutions

The equations of motion for the Linear and Quadratic Models are

Both models use the same Initial Conditions

DSolve gives explicit formulae for the velocity and position with the linear model.

The -> notation is a Mathematica replacement rule which facilitates substituting the solutions of this equation into other expressions.  The solutions can be readily extracted from the substitution rules "->" using the "/." command.  One slightly irritating thing is that due to the possibility of more than one solution and the extremely consistent output formats of Mathematica the solution is actually presented as a list of substitution rules.  Since we want the first (and only) solution we specify this with the Part command [[1]].

DSolve solves the Quadratic model but for some unobvious reason is unable to evaluate the constants corresponding to the intial conditions.  

Extracting the solutions from the substitution rules "->" using the "/." command again.

Note, the constants of integration (called C[1] and C[2]) in the solutions.  We can now evaluate at t = 0 and solve for C[1] and C[2] from the Initial conditions.

Which gives C[2] =0 and C[1] =  x0 and we can redefine the solutions with these values as follows.

The terminal velocities can be extracted directly from the equations (taking limits of the analytic solutions as t->∞ is possible - see below -  but in a more complicated situation might require specifying assumptions concerning the signs of g, k, m etc.) by setting v'[t] = 0 and solving for the velocity

Clearly in the quadratic cases we want the negative solution.

T-10 Model

Parameters

We first specify all the parameters required for the model.  Details are in the paper.

The terminal velocity of a sky-diver is approximately 100 miles per hour

Computation of Drag Coefficients

If the free-fall terminal velocity is known then setting the RHS of the equations of motion to zero and solving gives the drag coefficient.  The rate of descent for the T-10 parachute is 16 ft/se c[ebk, p.86] which gives the drag coefficient for the fully deployed parachute.  These could be defined as a function of the parameters above.

Analysis of the Model

Drag Coefficient

The quadratic model is complete specified once the drag coefficient is computed and the smoothness of the solution can be determined.  The following cell defines the drag coeeficient kk as a function of time "t" and substitutes in the appropriate values of the parameters defined above.

Graphs of kk and its derivative provide a first test of the smoothness of this function. However, be very careful before using a picture as a proof of continuity or smoothness.  Note, the drag coefficient during final descent is approximately 29. This value is a little lower than the prediction based on the terminal velocity.  This suggests that the landing velocity is likely to be a little higher than that for a jump from a 5-foot high wall and probably noticeably higher than the rate of descent published in[ebk,p.86].  The precise amount of this over-estimate remains to be seen.The drag coefficient during free-fall is much closer to the value obtained from the literature.

Note, the discontiunuity in the time rate of change of the drag coefficient.

A more reliable method of determining continuity of a function at a point is to compare its right-and left-hand limits at the point, i.e., the jump of the function at the point.   The drag coeeficient is continuous but not differentiable.  This means the acceleration will be continuous, but not the jerk.

Terminal Velocity

We are most interested in the free- fall and landing terminal velocities but the terminal velocity is easily found for any time during the jump.

The terminal velocity during free-fall and final descent are consistent with our earlier observations.The free-fall terminal velocity is extremely close to 100 mph.The final descent terminal velocity does not match reality nearlyy as well; as expected,the landing speed is higher than the speeds found in the literature.  Observe that the continuity of the drag coefficient says that the terminal velocity must be continuous for all time.This is not immediately apparent from the explicit definition of the terminal velocity or its graph.

Graphical Verification of Smoothness and Deployment

It is now easy to plug the explicit expression kk[t] in for k in hte quadratic model and substitute fro any remining paramters in the substitution list Param to get a readily solvable system of ODE.

Note the solution is returned as a list of InterpolatingFunctions as described earlier.  We exctract and plot the solutions.

Note the discontinuities in the jerk.  Now for a close-up.

Note the discontinuities in the jerk again.

Estimating the Landing Time

A complete plot of the descent is

The landing is near t=160 and we can obtain a better estimate using the built in Mathematica command FindRoot (which implements a variety of Newton like methods)

Apparently the landing occurs after 161.7 Seconds.  Checking this result