ORDINARY DIFFERENTIAL EQUATIONS POWERTOOLLesson 20 -- Application: Parachute ProblemProf. Douglas B. MeadeIndustrial Mathematics InstituteDepartment of MathematicsUniversity of South CarolinaColumbia, SC 29208
URL: http://www.math.sc.edu/~meade/E-mail: meade@math.sc.eduCopyright \251 2001 by Douglas B. MeadeAll rights reserved-------------------------------------------------------------------<Text-field layout="_pstyle14" style="_cstyle10">Outline of Lesson 20</Text-field>20.A The Parachute Problem20.A-1 Equations of Motion20.A-2 Terminal Velocity and Coefficients of Air Resistance20.A-3 The Full Model20.A-4 Explicit Solution20.A-5 Explicit Solution -- Simpler Form20.A-6 Numeric and Graphical Solution20.A-7 Time of Impact20.A-8 Final Comments20.A-9 References<Text-field layout="_pstyle14" style="_cstyle10">Initialization</Text-field>restart;with( DEtools ):with( plots ):<Text-field bookmark="20.A" layout="_pstyle14" style="_cstyle10">20.A The Parachute Problem</Text-field>Problem DescriptionTraining jumps for the Parachute Team at the United States Air Force Academy begin 4000 ft above ground level (AGL). The free-fall portion of the jump lasts about 10 seconds; free-falls longer than 13 seconds are grounds for removal from the team. The terminal velocity for free-fall is 120 mph. The landing velocity should be no worse than a free-fall from a 5' wall.How long after leaving the plane does a 190 pound skydiver return to solid ground?<Text-field bookmark="20.A-1" layout="_pstyle18" style="_cstyle15">20.A-1 Equations of Motion</Text-field>Under the assumption that the force of air resistance is proportional to velocity and opposes the motion, the second-order equation of motion for a parachutist falling in a coordinate system where NiMtJSJ4RzYjJSJ0Rw== is measured positive upward from the ground, iseom2 := m*diff( x(t), t$2 ) = -m*g - k*diff( x(t), t );with initial conditionsic2 := { x(0)=x0, D(x)(0)=0 };The equivalent first-order system of ODEs formed by writingNiMvLSUidkc2IyUidEcqJiUjZHhHIiIiJSNkdEchIiI=so that NiMvKiglImRHIiIjJSJ4RyIiIiokJSNkdEdGJiEiIiomJSNkdkdGKEYqRis=isodeX := diff( x(t), t ) = v(t);odeV := m*diff( v(t), t ) = -m*g - k*v(t);eomXV := { odeX, odeV }:The initial conditions for a parachutist stepping out of a plane at height NiMmJSJ4RzYjIiIh would beicX := x(0)=x0;icV := v(0)=0;icXV := { icX, icV }:<Text-field bookmark="20.A-2" layout="_pstyle18" style="_cstyle15">20.A-2 Terminal Velocity and Coefficients of Air Resistance</Text-field>The terminal velocity, NiMmJSJ2RzYjJSJURw==, is the equilibrium solution for the velocity ODE. That it, NiMmJSJ2RzYjJSJURw== must satisfyVterm := eval( odeV, v(t)=v[T] );because NiMqJiUjZHZHIiIiJSNkdEchIiI= approaches zero as NiMlInRH increases.The coefficient of air resistance is piecewise-defined with general formkk := piecewise( t<10, kFF, kL );The values for the drag coefficients during free-fall, NiMmJSJrRzYjJSNGRkc=, and landing, NiMmJSJrRzYjJSJMRw==, with units kg/s, can be determined from the information given in the problem. The values of the system parameters NiMlImdH, the gravitational constant, and NiMlIm1H, the mass of the parachutist, areparam0 := g = 9.81, m = 190 / 2.2;A free-fall velocity of 120 miles per hour, expressed in meters per second, isvFF := -120 * 5280*12*2.54/100/60/60;so the coefficient of air resistance during free-fall is found to bekFF := solve(eval( Vterm, [param0,v[T]=vFF] ), k );The value of NiMmJSJrRzYjJSJMRw== is based on the landing velocity for a jump from a 5' wall. Five feet becomex5 := 5 * 12*2.54/100;meters. The position and velocity for a free-fall jump from an initial height of NiMvJiUieEc2IyIiIS0lJkZsb2F0RzYkIiVDOiEiJA== m (= 5') aresol5 := eval( dsolve( eomXV union icXV, {x(t),v(t)} ), {param0, x0=x5, k=kFF} );The landing time, in seconds, for this short jump ist5 := fsolve( eval( x(t)=0, sol5 ), t=0..1 );The corresponding landing velocityv5 := eval( rhs( op(select(has,sol5,v(t))) ), t=t5 );is the final ingredient needed to determine the coefficient of air resistance during the time the parachute is deployedkL := solve(eval( Vterm, [param0,v[T]=v5] ), k );<Text-field bookmark="20.A-3" layout="_pstyle18" style="_cstyle15">20.A-3 The Full Model</Text-field>The complete model consists of the IVPsys := eval( eomXV, [param0, k=kk] );ic := eval( icXV, x0 = 4000 * 12*2.54/100 );Prior to obtaining an explicit solution to this IVP, note that the velocity satisfies a linear first-order ODE. The graph of the direction field for the velocity equation and the solution curve with NiMvLSUidkc2IyIiIUYn provides, in Figure 20.1, a first look at the solution to this model.DEplot( select(has,sys,diff(v(t),t)), {v(t)}, t=0..20, [[0,0]], v=-100..0, stepsize=0.01, linecolor=blue, title="Figure 20.1" );Observe how the velocity approaches the terminal free-fall velocity until the parachute is deployed at NiMvJSJ0RyIjNQ==. Thereafter, the velocity quickly approaches the new (and slower) terminal velocity.<Text-field bookmark="20.A-4" layout="_pstyle18" style="_cstyle15">20.A-4 Explicit Solution</Text-field>The explicit solution to this model issol := dsolve( sys union ic, { x(t), v(t) } );This solution is quite complicated even with floating-point coefficients, as obtained withevalf( sol );<Text-field bookmark="20.A-5" layout="_pstyle18" style="_cstyle15">20.A-5 Explicit Solution -- Simpler Form</Text-field>A slightly simpler representation of the solution can be obtained by solving the system twice, first with NiMvJSJrRyZGJDYjJSNGRkc= and the original initial conditions and then with NiMvJSJrRyZGJDYjJSJMRw== using the position and height from the free-fall solution at the time the parachute is deployed as initial conditions. The system with NiMvJSJrRyZGJDYjJSNGRkc= issys1 := eval( eomXV, [param0, k=kFF] );The explicit solution during free-fall issol1 := dsolve( sys1 union ic, { x(t), v(t) } );The equations of motion for the skydiver after the parachute is deployed aresys2 := eval( eomXV, [param0, k=kL] );The skydiver's height and velocity at the time the parachute is deployed areic2 := eval( {x(10)=x(10.),v(10)=v(10.)}, eval(sol1,t=10.) );The explicit solution after the parachute is deployed issol2 := dsolve( sys2 union ic2, { x(t), v(t) } ); Assembling the solution for both stages of the jump into a single expression leads tosolX := x(t) = piecewise( t<10, eval( x(t), sol1 ), eval( x(t), sol2 ) );for the position, andsolV := v(t) = piecewise( t<10, eval( v(t), sol1 ), eval( v(t), sol2 ) );for the velocity. As expected, this solution is in a much simpler form than the original solution.<Text-field layout="_pstyle21" style="_cstyle18">Note</Text-field>The explicit solution can be obtained in the same manner without specifying values for the parameters in the problem. This purely symbolic solution is quite complicated, but has a number of uses, including a sensitivity analysis of the solution with respect to each of the parameters.<Text-field bookmark="20.A-6" layout="_pstyle18" style="_cstyle15">20.A-6 Numeric and Graphical Solution</Text-field>A numeric solution is much easier to obtain and can be used to answer many questions about the motion. To confirm that Maple is computing a reasonable solution, recreate, in Figure 20.2, a plot of the velocity function over the first 20 seconds of the jump. (Since it is so easy to do, include a scaled version of the height.)solN := dsolve( sys union ic, { x(t), v(t) }, type=numeric ):odeplot( solN, [ [t,x(t)/10], [t,v(t)] ], 0..20, numpoints=500, legend=["x/10","v"], title="Figure 20.2" );<Text-field bookmark="20.A-7" layout="_pstyle18" style="_cstyle15">20.A-7 Time of Impact</Text-field><Text-field layout="_pstyle21" style="_cstyle18">Graphical Solution</Text-field>To see the motion for the entire jump, extend the time interval to 200 seconds as shown in Figure 20.3.odeplot( solN, [ [t,x(t)/10], [t,v(t)] ], 0..200, numpoints=500, legend=["x/10","v"], title="Figure 20.2" );This suggests that the skydiver lands in just over 3 minutes (180 seconds), as verified bysolN(180);A quick search leads to a better estimate of the impact time, namely,solN(181.76);Note that the impact velocity is essentially the same as the velocity for a free-fall from a 5' wall, which wasv5;<Text-field layout="_pstyle21" style="_cstyle18">Analytic Solution</Text-field>To determine the time of impact from the explicit solution, one might try solving viaq1 := solve( rhs(solX)=0, t );or in floating-point form,evalf( q1 );but this answer is not physically realistic. Further thought leads to the realization that this "solution" is obtained from the free-fall phase of the solution but the time of impact must occur during the post-deployment phase. To convey this additional knowlege to Maple, useTimpact := solve( eval( x(t), sol2 )=0, t );or in floating-point form,evalf( Timpact );and note the agreement with the time of impact obtained from the graphical and numerical solutions.<Text-field bookmark="20.A-8" layout="_pstyle18" style="_cstyle15">20.A-8 Discussion of the Model</Text-field>This version of the parachute problem is the classical model found in many ODE textbooks. It's major flaw is that the acceleration is discontinuous at the time the parachute is deployed.The acceleration can be obtained by differentiation of the velocity, givingsolA := eval( a(t) = diff( v(t), t ), diff( solV, t ) );The acceleration appears to have a jump discontinuity at the instant the ripcord is pulled, which is found to bet_jump := op( discont( rhs(solA), t ) );The size of the jump in the discontinuity isjump := evalf( limit( rhs(solA), t=10, right ) - limit( rhs(solA), t=10, left ) );`` = eval( jump/g, [param0] ) * g;An 8g acceleration is strong enough to (at least) seriously injure the skydiver. The physics of skydiving require the acceleration to be a NiMpJSJDRyUiMUc= function, i.e., the derivative of the acceleration, the jerk (think about it!), must be continuous. In the model analyzed in this worksheet, the acceleration is discontinuous! To obtain a physically realistic model it is necessary to model the (short) time during which the parachute is deployed. Models with this property are discussed in the reference articles authored by D. B. Meade and A. A. Struthers.The references discuss several interesting problems based on the basic parachute problem. For example, the article by J. Drucker and the D. B. Meade's MapleTech article discuss finding the deployment time for the shortest jump that lands with an impact velocity within, say, 10% of the typical impact velocity.<Text-field bookmark="20.A-9" layout="_pstyle18" style="_cstyle15">20.A-9 References</Text-field>J. Drucker, Minimal time of descent, The College Math. J., 26 (1995), pp. 232-235.D. B. Meade, Maple and the parachute problem: modelling with an impact, MapleTech, 4 (1997), pp. 68-76.D. B. Meade, ODE models for the parachute problem, SIAM Review, 40 (1998), pp. 252-255.D. B. Meade and A. A. Struthers, Differential equations in the new millenium: the parachute problem, Int. J. Engng. Ed., 15 (1999), pp. 417-424.Note:For Maple worksheets and other materials related to the papers (co-)authored by D. B. Meade, please see his list of publications at http://www.math.sc.edu/~meade/publ.html.[Back to ODE Powertool Table of Contents]