REVITALIZING DIFFERENTIAL EQUATIONS WITH A COORDINATED USE OF NUMERICAL, GRAPHICAL, AND ANALYTICAL METHODS

Douglas B. Meade Department of Mathematics University of South Carolina Columbia, SC 29208

                          meade@math.sc.edu

Introduction

The typical introductory course in differential equations often degrades into the discussion of a collection of solution techniques. While a thorough knowledge of these techniques is useful, the course should also develop an appreciation for the general theory and application of differential equations. Textbooks and course materials that attempt to emphasize this broader approach to differential equations have appeared in recent years ([1], [2]). The same effect can be achieved without changing the textbook; all that is required is a little creativity from the instructor (and the students).

The objective of this paper is to demonstrate how a traditional problem from a traditional text can be modified to provide a more realistic view of the applicability of differential equations. Almost any “applied” problem – compartment models, mechanics, population dynamics, Newton’s law of cooling, electrical circuits – from a traditional text could be selected as the basis for this discussion. The specific problem selected for this paper is based on the application of Newtonian mechanics to a body falling under the influence of gravity and air resistance [3, p. 94].

A parachutist whose mass is 75 kg drops from a helicopter hovering 4000 m above the ground and falls towards the earth under the influence of gravity. Assume that the gravitational force is constant and that the force due to air resistance is proportional to the velocity of the parachutist, with the proportionality constant k = 15 kg/sec when the chute is closed and with constant k = 105 kg/sec when the chute is open. If the chute does not open until 1 min after the parachutist leaves the helicopter, after how many seconds will she hit the ground?

The piecewise-defined resistance coefficient presents an interesting twist to this problem. This problem, which initially appears to be a simple second-order initial value problem (IVP), is easier to solve as a system of two first-order IVPs. (Note that the discontinuity occurs in the coefficient of the velocity, not in the forcing term. Thus, Laplace transforms are not as natural a choice as might be expected.)

Many extensions of this problem are possible. One family of variations can be created by changing the criteria for the opening of the chute. For example, instead of deploying after one minute, answer the same questions when the chute opens at a given altitude or when the parachutist reaches a prescribed velocity. The extensions considered in this paper are based around the impact velocity. To be more specific, let t_{1} denote the time when the chute opens and k_{1} and k_{2} denote the resistance coefficients during free-fall and after the chute has deployed, respectively. The specific enhancements discussed in this paper are:

• answer the same question when the resistance coefficient after the chute opens is k_{2} = 30 kg/sec (instead of 105 kg/sec). Compare the impact velocities.

• for what values of k_{2} will the impact velocity be no greater than 10 m/sec?

• what is the longest free-fall that can be made that will ensure that the impact velocity is no more than 7.25 m/sec (use k_{1} = 15 m/sec and k_{2} = 105 m/sec).

Each of these problems is increasingly more difficult. The solutions will be obtained using a combination of analytic, graphical, and numerical routines available in Maple. The first step is to answer the original question.

The Solution of the Original Problem

This situation can be modelled either in terms of a second-order differential equation in terms of the distance fallen, x,

 > ODE := m*diff( x(t), t$2 ) = m*g - k*diff( x(t), t );

ODE := m (\frac{∂^{2}}{∂t^{2}} x( t)) = m g − k (\frac{∂}{∂t} x( t))

or as a system of two first-order equations involving x and v, the velocity (down is the positive direction),

 > SYS := [ diff( x(t), t ) = v(t), m*diff( v(t), t ) = m*g-k*v(t) ];

SYS := [\frac{∂}{∂t} x( t) = v( t), m (\frac{∂}{∂t} v( t)) = m g − k v( t)]

For the original problem the coefficients are

 > PARAM := m=75,k=k1+(k2-k1)*Heaviside(t-t1),k1=15,k2=105,t1=60,g=9.81:

A graphical solution (see Figure 1), provides a good first look at the solution.

 > with( DEtools ):
 > DEplot( subs( PARAM, SYS[2] ), [t,v], 0..250, { [0,0] }, v=0..60,
 >  stepsize=1, title=`Direction Field: velocity vs time` );
 > DEplot( subs( PARAM, SYS ), [t,x,v], 0..250, { [0,0,0] }, x=0..4000,
 >  scene=[t,x], stepsize=1, title=`Solution: distance vs. time` );

Rough estimates of the answers to the original problem can be deduced from these plots. The parachutist returns to earth in approximately 240 sec. The velocity at impact is approximately 7 m/sec. When the parachute opened, after 1 min, the velocity was approximately 49 m/sec and the parachutist had fallen about 2700 m.

Figure 1: The velocity and distance fallen for k_{1} = 15 kg/sec and k_{2} = 105 kg/sec.

The explicit solution is not difficult to find, and will yield more accurate answers. The differential equation that models the first minute of the jump is

 > SYS1 := subs( k=k1, {op(SYS)} );

SYS1 := \{\frac{∂}{∂t} x( t) = v( t), m (\frac{∂}{∂t} v( t)) = m g − k1 v( t)\}

with initial conditions

 > IC1 := { x(0)=0, v(0)=0 };

IC1 := { v( 0 ) = 0, x( 0 ) = 0 }

This initial value problem is easily solved by hand. The velocity equation is firstorder, linear and independent of x; solve for v, then integrate to find x.

 > SOL1 := dsolve( SYS1 union IC1, {x(t),v(t)} ):
 > X1 := subs( SOL1, x(t) ): V1 := subs( SOL1, v(t) ):
 > 'X1' = X1, 'V1' = V1;

X1 = \frac{m g t k1 − m^{2} g + m^{2} g e^{(− \frac{k1 t}{m})}}{k1^{2}} , V1 = \frac{m g}{k1} − \frac{m g e^{(− \frac{k1 t}{m})}}{k1}

The second stage of the jump is similar,

 > SYS2 := subs( k=k2, {op(SYS)} );

SYS2 := \{m (\frac{∂}{∂t} v( t)) = m g − k2 v( t), \frac{∂}{∂t} x( t) = v( t)\}

except that the “initial conditions” must ensure continuity with the free-fall solution

 > IC2 := { x(t1) = subs(t=t1,X1), v(t1)=subs(t=t1,V1) };

IC2 := \{x( t1 ) = \frac{m g t1 k1 − m^{2} g + m^{2} g e^{(− \frac{k1 t1}{m})}}{k1^{2}} , v( t1 ) = \frac{m g}{k1} − \frac{m g e^{(− \frac{k1 t1}{m})}}{k1}\}

The solution is obtained exactly as before, except for algebraic complications.

V2 := \frac{m g}{k2} − \frac{m g (k1 − k2 + e^{(− \frac{k1 t1}{m})} k2 ) e^{(− \frac{k2 t}{m})}}{e^{(− \frac{k2 t1}{m})} k2 k1}

> SOL2 := dsolve( SYS2 union IC2, {x(t),v(t)} ):

> X2 := subs( SOL2, x(t) ): V2 := subs( SOL2, v(t) );

The complete solution can now be expressed in terms of the Heaviside function

 > X := X1 + (X2-X1)*Heaviside(t-t1):
 > V := V1 + (V2-V1)*Heaviside(t-t1):

The total time of the jump and the landing velocity are now easily determined. The exact values can be determined, but only in terms of Lambert’s W function.

 > t_impact := fsolve( subs( PARAM, X2 )=4000, t );
 > v_impact := evalf(subs( t=t_impact, PARAM, V ) );

t impact := 241.5601768
v impact := 7.00714285

The Three Enhancements

The first problem can be answered by simply changing the parameter values

 > PARAM := m=75,k=k1+(k2-k1)*Heaviside(t-t1),k1=15,k2=30,t1=60,g=9.81:

and re-computing the time and velocity at impact

 > t_impact := fsolve( subs( PARAM, X2 )=4000, t );
 > v_impact := evalf(subs( t=t_impact, PARAM, V ) );

t impact := 110.5988480
v impact := 24.52500004

The second problem requires a little more thought. The fact that the time of impact depends on the value of k_{2} makes it difficult to determine when the impact velocity will not exceed a prescribed threshold, v_{max}.. Here is where the graphical solutions are useful: the velocity decreases exponentially to its terminal velocity v_{∞} = \frac{mg}{k_{2}}, once the chute is deployed. Thus, unless the chute opens immediately prior to impact, the landing and terminal velocities are nearly equal. Choosing k_{2} > \frac{mg}{v_{max}} ensures that v_{∞} < v_{max}. Since the parachutist is still more than 1000 m above ground when the chute opens, the above condition on k_{2} also ensures that the impact velocity will not exceed v_{max}.

And, finally, the control problem. The objective is to find the latest time, t_{1}, at which the parachute can be opened and still ensure an impact that does not exceed 7.25 m/sec. Using the original values of all parameters (except t_{1}), the velocity and distance fallen after the chute opens are given by (output suppressed):

 > PARAM := m=75, g=9.81, k1=15, k2=105:
 > XX := fnormal( subs( PARAM, X2 ), 4 ):
 > VV := fnormal( subs( PARAM, V2 ), 4 ):

In general, many pairs (t_{1}, t) satisfy XX(t; t_{1}) = 4000 and V V (t; t_{1}) ≤ v_{max}. Here,

 > Vmax := 7.25:

In general, for each value of t there will be a different latest time, t_{1}, at which the chute can be opened without violating the prescribed conditions. The most extreme case would be to take t_{1} ≥ t (i.e., never open the chute). The time, t_{ff}, required to free-fall 4000 m, and the corresponding impact velocity, are found to be

 > t_ff := fsolve( subs( PARAM, X1 )=4000, t=0..300 );
 > v_ff := evalf( subs( PARAM, t=t_ff, V1 ) );

t ff := 86.54943920

v ff := 49.04999851

It is not surprising that the landing velocity exceeds the threshold. This test does, however, yield enough information to narrow the search to t_{1} < t_{ff} < t. To illustrate, consider the two cases when impact occurs after 88 and after 100 seconds.

 > plot( subs( t=88, { XX-4000, VV-Vmax } ), t1=80..t_ff,
 >  -5..5, title=`Case 1: t=88 seconds` );
 > plot( subs( t=100, { XX-4000, VV-Vmax } ), t1=80..t_ff,
 >  -5..5, title=`Case 2: t=100 seconds` );

Figure 2: Plots of XX(t; t_{1})−4000 and V V (t, t_{1})−7.25 for t = 88 and t = 100 seconds.

This plots show that an impact time of 88 seconds (t = 88) is achieved when t_{1} ≈ 85.6 sec; the corresponding impact velocity exceeds the prescribed threshold by more than 1 m/sec. The 100 second jump calls for chute deployment only two seconds earlier, t_{1} ≈ 83.6 sec, but the landing occurs at a velocity safely under the threshold. The shortest jump which meets the constraints can be found through additional tests.

 > Tset := [ 86.6, 87, 88, 89, 90, 100 ]:
 > T1set := [ seq( [ T, fsolve( subs(t=T,XX=4000), t1 ) ], T=Tset ) ]:
 > Vset := [ seq( fnormal(subs(t=p[1],t1=p[2],[t,t1,VV]), 4), p=T1set ) ];

Vset := [[ 86.6, 86.31, 35.19 ], [ 87., 85.92, 16.39 ], [ 88., 85.62, 8.558 ],
[ 89., 85.44, 7.290 ], [ 90., 85.27, 7.064 ], [ 100., 83.60, 7.007 ]]

Note that the shortest possible jump from 4000 m with an impact velocity not exceeding 7.25 m/sec would last slightly more than 89 seconds.

References

[1] Paul Blanchard, Robert Devaney, Glen Hall, Differential Equations, PWS Publishing Co., 1996.

[2] Robert L. Borrelli, Courtney Coleman, William E. Boyce, Differential Equations Laboratory Workbook: A Collection of Experiments, Explorations and Modeling Projects for the Computer, John Wiley & Sons, Inc., 1992.

[3] R. Kent Nagel and Edward B. Saff, Fundamentals of Differential Equations, 3rd ed., Addison-Wesley, 1993.