Parachute Problem

provide a good description of the motion of the skydiver. Given an initial value problem in terms of x and v, the procedure motion uses dsolve,numeric to compute the motion of the object at user-specified times; plot_motion uses the information computed by motion to creates a single plot containing any combination of the position, velocity, acceleration, and jerk.

> motion := proc( IVP:set(equation),
>  TIMES:list(numeric), t )
> local POS, VEL, ACC, JERK, pts, acc, jerk,
>  soln, T;
> options `Copyright 1996 by Douglas B. Meade`;
> soln:=dsolve(IVP, {x(t),v(t)}, type=numeric,
>  output=listprocedure):
> POS:=subs( soln, x(t) );
> VEL:=subs( soln, v(t) );
> acc:=solve( ode(IVP,v), diff(v(t),t) );
> ACC:=unapply( subs( v(t)='VEL(t)', acc ), t );
> jerk:=diff(acc,t);
> JERK:=unapply(subs( {v(t)='VEL(t)',
>  'diff(v(t),t)'='ACC(t)'},
>  jerk ), t );
> pts:=seq(
>  evalf([T,POS(T),VEL(T),ACC(T),JERK(T)]),
>  T=sort(TIMES) ):
> RETURN( [ pts ] );
> end:
> plot_motion :=
> proc( pts:list, scale:set(equation) )
> local P, PLOTS, S, VARS;
> options `Copyright 1996 by Douglas B. Meade`;
> PLOTS := NULL;
> VARS := map(lhs,scale);
> if member(X,VARS) then
>  S:=subs(scale,X);
>  PLOTS:=PLOTS,
>  plot([ seq( [P[1],P[2]/S], P=pts )],
>  color = RED ):
> fi;
> if member(V,VARS) then
>  S:=subs(scale,V);
>  PLOTS:=PLOTS,
>  plot([ seq( [P[1],P[3]/S], P=pts )],
>  color = BLUE ):
> fi;
> if member(A,VARS) then
>  S:=subs(scale,A);
>  PLOTS:=PLOTS,
>  plot([ seq( [P[1],P[4]/S], P=pts )],
>  color = GREEN ):
> fi;
> if member(J,VARS) then
>  S:=subs(scale,J);
>  PLOTS:=PLOTS,
>  plot([ seq( [P[1],P[5]/S], P=pts )],
>  color = CYAN ):
> fi;
> RETURN( display( { PLOTS } ) );
> end:

Using Maple’s Heaviside function to implement the piecewise-constant air resistance (2),

> K[pw_const] :=  k1+(k2-k1)*Heaviside(t-td):

the motion of the skydiver is governed by the initial value problem

> IVP := subs( k=K[pw_const], SYS ) union IC:

The specific motion can be produced after specifying values for each parameter. The values given in [9] are

> PARAM[ns] := k1=15, k2=105, td=60, x0=4000,
>  m=75, g=981/100:

Lastly, the list of times is selected with a finer resolution near the time of deployment:

                                                        > snapshot[ns] := [ $(0..58),
                                                        >  i/10 $ i = 590 .. 649,
                                                        >  5*i $ i = 13 .. 60 ]:

An informative graph of the motion during the first five minutes of the jump is obtained as follows:

                                                        > PTS[ns]:=motion( subs(PARAM[ns],IVP),
                                                        >  snapshot[ns], t ):
                                                        > scale[noJ] := { 'X'=1000, 'V'=10, 'A'=981/100 }:
                                                        > display(
                                                        > {
                                                        >  plot_motion( PTS[ns], scale[noJ] ),
                                                        >  textplot( [ [ 120, 1.5, `x/1000` ],
                                                        >  [ 80, -1, `v/10` ],
                                                        >  [ 70, 5, `a/g` ] ] )
                                                        > },
                                                        > title=`Motion for Nagle/Saff parameters` );

This plot supplies good first approximations to a number of interesting characteristics of the jump. When the ripcord is pulled at the end of the first minute, the skydiver is still 1300 m above ground level and descending at an almost constant velocity of 49 m/s. Within 5 seconds from the time the parachute opens the motion settles into an essentially steady descent at 7 m/s; the total length of the jump is very close to 4 minutes.

More precise answers can be obtained from an explicit solution to (1). The solution during the free-fall portion of the jump is:

                                                        > SOL1 := dsolve( subs(k=k1,SYS) union IC,
                                                        >  { x(t), v(t) } ):
                                                        > V1 := subs( SOL1, v(t) );
                                                        > X1 := subs( SOL1, x(t) );

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

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

/k1^{2}

“Initial conditions” for the second stage of the jump are

\overline{2}