The Units Package with Derivatives, Integrals and DEs

The new Units package in Maple 7 can incorporate units in solutions to problems involving addition, multiplication and division of physical quantities. A demonstration of this basic functionality is available at http://www.adeptscience. co.uk/go?pg=H34 However, the Units package also keeps track of units throughout computations involving derivatives, integrals and differential equations. In this article, we show you how.

In the examples below, we use the Standard Units environment, which allows you to name quantities without worrying about overwriting common symbols used to name units, such as m for "metre" and s for "second". For example, you could specify a speed of 5 metres per second with s:= 5*Unit(m/s), with no conflict involving the symbol s. The alternative to the Standard Units environment is the Natural Units environment, which is used in the Web demo referenced above.
> with( Units[Standard] ):

Example 1: Units with Calculus
From a tower 20 metres above the ground, you throw a ball at speed of 15 miles per hour, at an angle of 65¾ with the horizontal. How fast is the ball falling when it hits the ground? Give the answer in metres per second and in miles per hour.

We first enter the four known parameters in their respective units:
(1) The initial height. We use the syntax Unit(m) to specify that the height is in metres. The Maple output reflects the unit.
> h0 := 20 * Unit( m );

(2) The initial speed in miles per hour. Notice that Maple converts the speed to m/s by default. (The default unit system is SI but can be set by the user to any of eight systems recognized by the Units package, including FPS, CGS and atomic.)
> speed := 15 * Unit( mi/hour );


(3) The angle in degrees.
> theta := 65.0 * Unit(degrees);


(4) The constant acceleration of gravity in m/s2.
> a := -9.81 * Unit( m/s^2 );


Next, we compute the initial vertical velocity by multiplying the speed by sin(q). The Units package overrides Maple's assumption that arguments to sin() are in radians, allowing us to enter an angle in any angular unit, in this case, degrees.
> v0_up := speed * sin( theta );


And now for a little calculus. To get the ball’s falling speed on impact, we need a formula for the height as a function of time. We obtain this by two integrations. We obtain the vertical-velocity formula by integrating the acceleration with respect to time t. To integrate with respect to a variable having a unit — in this case, seconds – use the syntax int( integrand, variable * Unit(unit)). We supply the initial velocity as the integration constant. The output is given as a function of t in m/s.
> v := v0_up+int(a, t*Unit(s));

We can now compute the height formula as a function of time by integrating the velocity in t and adding the height of the tower as the integration constant. The output is given as a function of t in metres.
> h := h0+int( v, t*Unit(s));

To compute the time of impact, we solve h(t) = 0 for t.
> fallTime := solve( h = 0, t );


Finally, we can answer the question by evaluating the vertical velocity at the positive solution listed above using the eval command.
> eval( v, t=fallTime[1] );


The answer above is in metres per second, which we convert to mph using the convert command.
> convert( %, units, mi/hour );

Example 2: Units with Differential Equations
A mass of 23.1 grams is suspended from the ceiling by a spring of stiffness 0.1 kg/s2. Starting at 5.25 cm above its equilibrium, the mass is pushed upwards with an initial speed of 1 m/s. What will the displacement of the spring be in 10 seconds?
We first write down the DE governing the motion of an undamped spring-mass system with mass m and spring stiffness k.
> springDE := -k*x(t) = m*diff(x(t),t,t);


We specify the initial displacement and velocity as x0 and v0, respectively.
> springIC := x(0)=x0, D(x)(0)=v0;

Now we solve the DE with the dsolve command to obtain the displacement of the mass as a function of t in terms of the four parameters.
> displacement := dsolve( {springDE, springIC},x(t));

To incorporate units into the solution, we specify m, k, v0 and x0 in a list in their respective units. We also specify that t is measured in seconds, using the syntax t = t * Unit(s).
> springValues := [m = 23.1 * Unit( g ),
k = 0.1 * Unit( kg/s^2 ),
v0 = 1.0 * Unit( m/s ),
x0 = 5.25 * Unit( cm ),
t = t * Unit( s )];

Compute the displacement as a function of t alone by evaluating it on those parameters using the eval command. Notice that Maple has cancelled some units out of the numerators and denominators of the displacement formula, leaving length as the remaining dimension, as expected.
> eval( displacement, springValues ) ;

Finally, answer the question by evaluating the above at t=10. We don't need to enter t=10*Unit(s) because we already specified the units of t in the parameter list.
> eval( %, [t = 10] ) ;

Example 3: Energy Conversions
In the physical sciences, one must often compute energy conversions. Tables of energy conversion factors are given in major books of tables, for example, the CRC Handbook of Chemistry and Physics, and online at sites like http://orgwww.chem.uva.nl/cgi-bin/conv/. Perhaps counterintuitively, quantities measured in units of length, mass, frequency and temperature can, in fact, be converted to quantities measured in units of energy. Three classic energy equations of quantum physics enable such conversions:
E = mc2
E = hn
E = kT

In the above, m stands for mass, c for the speed of light, h for Planck's constant, n for frequency, k for Boltzmann's constant and T for temperature.

The Units package supports these energy conversions through the energy option in the convert command. For example, to convert one Hertz, one kilogram and one degree Kelvin to Joules (the SI unit of energy), issue the following command:
> map(convert, [Unit(1/s), Unit(kg), Unit(K)], units, J, energy);

Using energy as the "medium of conversion", all the units mentioned above that are convertible to energy are also inter-convertible. For example, we can convert a mass of 0.732522e-37 kilograms to a wavelength in 1/centimetres as follows:
> convert( 0.732522e-37*Unit(kg), units, 1/cm, energy );

Thus, 0.732522e-37 kilograms is approximately equivalent to 331 wavenumbers.
Where Maple stands out from the crowd is in its ability to mix standard conversions and energy conversions correctly. For example, consider the following conversion of a force of 0.732522e-37 Newtons to Watts:
> convert( 0.732522e-37*Unit(N), units, W, energy );


To compute this using standard printed conversion tables, it would first be necessary to determine that Newtons (N) = Joules/meter and Power (P) = Joules/second. Thus, the above conversion is identical to converting wave numbers to frequency, as demonstrated by the following conversion of 0.732522e-37 wavenumbers to Hertz:
> convert(0.732522e-37*Unit(1/m),units,1/s, energy);