Application Center - Maplesoft

App Preview:

Units and dimensions in problems of calculus and ODEs

You can switch back to the summary page by clicking here.

Learn about Maple
Download Application


 

UnitsCalculus.mws

The Units Package with Derivatives, Integrals and ODEs

by Waterloo Maple

The 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.mapleapps.com/categories/maple7/html/Units.html.) However, the Units package also keeps track of units throughout computations involving derivatives , integrals and differential equations . In this application, 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 "meter" and s for "second". For example, you could specify a speed of 5 meters 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.

The Standard Units environment specifies that a speed of 5 meters per second would be entered in Maple as 5*Unit(m/s) , as opposed to 5*m/s , as would be done using the Natural Units environment.

> restart: with(Units[Standard]):

Units with Calculus: The Falling Body Problem

From a tower 20 meters above the ground, you throw a baseball at speed of 15 miles per hour, at an angle of 65 with the horizontal. How fast is the baseball falling when it hits the ground? Give the answer in meters 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 meters. The Maple output reflects the unit.

> h0 := 20 * Unit( m );

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 );

speed := 4191/625*Unit([m/s])

(3) The angle in degrees

> theta := 65.0 * Unit( degrees );

theta := 65.0*Unit([arcdeg])

(4) The constant acceleration of gravity in meters per second squared . We assume that the y-axis points upwards, so that the acceleration is negative.

> a := -9.81 * Unit(m/s^2);

a := -9.81*Unit([m/s^2])

Next, we compute the initial vertical velocity by multiplying the speed by sin(theta) . 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 );

v0_up := 6.079*Unit([m/s])

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));

v := (6.079-9.810*t)*Unit([m/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 meters.

> h := h0 + int(v, t*Unit(s));

h := (20+6.079*t-4.905*t^2)*Unit([m])

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

> fallTime := solve(h = 0, t);

fallTime := 2.732, -1.493

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] );

-20.72*Unit([m/s])

The answer above is in meters per second, which we convert to miles per hour using the convert command.

> convert(%, units, mi/hour);

-46.35*Unit([mi/h])

Units with ODEs: The Spring Problem

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 ODE 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);

springDE := -k*x(t) = m*diff(x(t),`$`(t,2))

We specify the initial displacement and velocity as x0 and v0 , respectively

> springIC := x(0)=x0, D(x)(0)=v0;

springIC := x(0) = x0, D(x)(0) = v0

Now we solve the ODE 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) );

displacement := x(t) = v0*m^(1/2)/k^(1/2)*sin(1/m^(1/2)*k^(1/2)*t)+x0*cos(1/m^(1/2)*k^(1/2)*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 )];

springValues := [m = 23.1*Unit([g]), k = .1*Unit([kg/s^2]), v0 = 1.0*Unit([m/s]), x0 = 5.25*Unit([cm]), t = t*Unit([s])]

> displacement, springValues;

x(t) = v0*m^(1/2)/k^(1/2)*sin(1/m^(1/2)*k^(1/2)*t)+x0*cos(1/m^(1/2)*k^(1/2)*t), [m = 23.1*Unit([g]), k = .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 canceled some units out of the numerators and denominators of the displacement formula, leaving length as the remaining dimension, as expected.

> eval( displacement, springValues ) ;

x(t*Unit([s])) = (.4806*sin(2.081*t)+.5250e-1*cos(2.081*t))*Unit([m])

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] ) ;

x(10*Unit([s])) = .4247*Unit([m])