Auxiliary Procedure Definitions (must be executed -- just press return)

As the system of ODEs and initial conditions are both sets, it will be useful to have a means to extract a specific entry from either object.

> ic := (IVP,y0) -> op( select( has, IVP, y0 ) ):
ode := (IVP,y) -> op( select( has, IVP, diff(y(t),t) ) ):

>

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. Given an initial value problem for the rate of change of the velocity and the initial position, the procedure motion uses dsolve,numeric to compute a numerical approximation to the position, velocity, acceleration, and jerk of the object. The result of a call to motion is a list of procedures. plot_motion uses the information computed by motion to create a single plot containing any combination of the position, velocity, acceleration, and jerk. The plot returned by plot_motion is created with odeplot .

Note that it may be necessary to exercise special care near discontinuities of the solution and/or ODE. In order to "see" discontinuities in the solutions, a point plot is recommended; to obtain a point plot, simply add style=POINT as an optional argument to the plot_motion command. For finer resolution, I also add numpoints=200 .

>

> motion := proc( IVP::set(equation), t::name )
local ODEpos, auxIVP, soln;
options `Copyright 1999 by Douglas B. Meade`;
ODEpos:=diff( x(t), t ) = v(t);
auxIVP := { ODEpos };
soln:=dsolve(IVP union auxIVP, {x(t),v(t)},
type=numeric, output=listprocedure ):
RETURN( soln );
end:

> plot_motion := proc( soln::list, time_int::{range(numeric),list(range(numeric))}, scale::set(equation) )
local acc, jerk, pos, vel, CURVES, S, TI, TIME_INT, opts;
options `Copyright 1999 by Douglas B. Meade`;
if nargs>3 then opts:=args[4..-1] else opts:=NULL fi;
if map(lhs,scale) minus {x,v,a,j} <> {} then
ERROR("invalid argument; scale factors can be specified only for x, v, a, and j, but received",scale);
fi;
soln(0); # always begin the plot by integrating from the initial point
pos := x=x(t);
vel := v=v(t);
acc := a=solve( ode(IVP,v), diff(v(t),t) );
jerk := j=rhs( diff( acc, t ) );
jerk := subs( 'diff(x(t),t)'=rhs(vel), 'diff(v(t),t)'=rhs(acc), jerk );
CURVES := subs( {pos, vel, acc, jerk}, [seq( [t,lhs(S)/rhs(S)], S=scale)] );

RETURN( plots[display]( [seq( plots[odeplot]( soln, CURVES, TI, opts ),
TI=convert(time_int,list) )] ) );
end:

>

The ``jump'' in a function at a given point is the difference between the left and right limits of the function at that point.

> jump := ( f, pt ) -> limit( f, pt, right ) - limit( f, pt, left ):

>