bump85.gif

A System with Air Resistance

The Fall of The Rubber Duck Revisited Tutorial Index
An Introduction to Differential Equations Maple Index
First-Order Homogeneous Linear Diff Eqs How to...
First-Order Nonhomogeneous Linear Diff Eqs Problem Set Index

You can stop holding your breath, because we're leaving fantasy land and entering the real world, where there is plenty of air to breath and complicate our systems.  At this point most students typically have enough physics knowledge to understand the basic concepts of these more realistic systems, but their mathematical toolboxes lack a few of the necessary tools for handling the associated mathematical models.  The missing tools are a knowledge of differential equations and the methods for solving them.  Fortunately, for the type of problem introduced in this section, the requisite knowledge is not difficult to learn.  If you are at least concurrently taking a course in integral calculus, and can integrate the function 1/x , then you know enough to comprehend and benefit from the material presented in this section.

Read the material linked to by the following three hyperlinks from top down.  Together they provide introductory information about differential equations and their solution techniques.  After reading these you will be well prepared to handle the discussion accompanying the example physics problem below.

An Introduction to Differential Equations
First-Order Homogeneous Linear Differential Equations
First-Order Nonhomogeneous Linear Differential Equations

New Maple Operators

  • dsolve   The ordinary differential equation solver.  Some examples that we solved by hand in the differential equation introduction:

    > diff_eq1 := diff(y(t),t)+2/t*y(t) = 0;

    diff_eq1 := diff(y(t),t) + (2/t)*y(t) = 0

    > soln1 := dsolve( diff_eq1, y(t) );

    soln1 := y(t) = (_C1)/t^2

    The variable _C1 is a Maple produced constant.  You can check that this is the same answer as we computed by hand.

    > diff_eq2 := diff(y(t),t)-4*y(t) = 7;

    diff_eq2 := diff(y(t),t) - 4 y(t) = 7

    > soln2 := dsolve( diff_eq2, y(t) );

    soln2 := (-7/4) + _C1 e^(4 t)

    Again, you can compare this solution to the one we computed by hand.

    For an example of how to set up dsolve for an initial-value problem, see the example physics problem below.  Also, check out ?dsolve.

  • odetest   This function checks to see that a given equation satisfies a specific ordinary differential equation.  A zero result is a confirmation, but a nonzero result does not necessarily mean that the equation is not a solution. ?odetest.

    > odetest( soln1, diff_eq1 );

    0

    > odetest( soln2, diff_eq2 );

    0

    You can see the by hand check of soln1 here, and likewise for soln2 here.

  • for  do   Maple is a powerful scripting language.  It contains all the basic control structures that you would expect of a full featured scripting language, including a for looping structure.

    > for i from 1 by 2 to 7 do i = i*i od;

    1 = 1
    3 = 9
    5 = 25
    7 = 49

    > for j from 1 to 4 do j = 2*j od;

    1 = 2
    2 = 4
    3 = 6
    4 = 8

    The variable following for is the loop index.  Without the by clause, the default increment value is one.  The statements surrounded by do and od are executed for each value of the index.  Also check out ?seq (sequence function), and for more on Maple programming start with ?statements.

Top

horzn_ln.gif (2407 bytes)

The Fall of the Rubber Duck Revisited

Let a rubber duck fall from the top of a 500 m cliff just as in The Fall of the Rubber Duck, but this time in the analysis of the situation, include the drag force induced on the rubber duck by its movement through the air.  This additional consideration results in some significant changes in the mathematical model of the problem.  In this analysis we will also need to know the rubber duck's mass; let it be 3 kg.

It turns out that the drag force on an object is usually directly proportional to its velocity, D = cV, where c is some constant, or the square of the velocity, D = cV^2.   This fact was determined experimentally, not by theoretical deduction or mathematical proof.  The value of c is also determined experimentally, and is affected by things such as an object's surface composition, length, general shape, etc.

The Reynolds number helps to determine which proportionality to use.  At any given velocity,

R = ( rho V L)/ mu

where L is the length of the object, rho the density, and mu the viscosity of the medium through which the object travels.  If R < 10 the proportionality is linear, if R > 1000, the proportionality is to V^2, and if 10 < R < 1000, neither linearity or the square reflects the situation.  Let's assume that the rubber duck is small and light enough that it reaches what is called its terminal velocity before the Reynolds number has the chance to grow larger than 10.  This will keep the differential equation that describes the situation linear (and thus we know how to solve it).

At a typical moment in the duck's fall, the force diagram might look like

duck force diagram

From the diagram and Newton's second law we see that at any given moment during the fall,

F[net]= D - W

While this equation may look a lot like many others you have seen, it's very likely that it is quite different in that a force on its right side is not a constant with respect to time; D varies with time.  Newton's second law and component forces that vary with time combine to create not your run of the mill algebraic equations, but differential equations, and you must use the solution methods of differential equations to extract useful information from them.

Note that,

F[net] = m diff( y(t), t, t)

D = cV

or,

D = c diff( y(t), t )

W = m g

Thus,

m diff( y(t), t, t ) = c diff( y(t), t ) - m g

Rearranging it into standard form,

diff( y(t), t, t ) - (c/m) diff( y(t), t) = -g

Of course, this is a second-order nonhomogeneous linear differential equation, and we've only studied how to solve first order equations of this type.  But, the following observations enable us to convert this into a first-order nonhomogeneous linear differential equation.

v(t) = diff( y(t), t )

And,

diff( v(t), t ) = diff( y(t), t, t )

Thus,

diff( v(t), t ) - (c/m) v(t) = -g

Top

horzn_ln.gif (2407 bytes)

At this point we begin a Maple solution to this initial-value problem.

> restart;

> diff_eq := diff(v(t),t)-v(t)*(c/m)=-g;

diff_eq := diff( v(t), t ) - (c/m) v(t) = -g

> dsolve({diff_eq, v(0)= 0}, v(t) );

v(t) = (g m/c) - (g m/c) e^(c t/m)

> odetest(%,%%);    # a check

0

> vel := simplify(%%);     Note the simplicity of this solution as opposed to the by hand method of the differential equation introduction; three in-depth pages of information reduced to just three commands.

vel := v(t) = - (g m/c) ( -1 + e^(c t/m) )

Keep in mind that c is not a constant produced during the solution process by any indefinite integrals, but the constant of proportionality from the equality D = cV (this is an initial-value problem, and thus includes only definite integrals in the solution).  In most situations, the magnitude of this constant is between 0.4 and 1.0.  Our choice of orientation for the vertical axis implies that our constant will need to be negative (see problem 1 of Differential Equation Problem Set).  Let c = -0.8, m = 3 kg , and g = 10 m / s ^2.

> c := -0.8; m := 3; g := 10;

c:= - .8
m := 3
g := 10

> vel;

v(t) = -37.50... + 37. 50... e^(-.266... t)

> Digits := 3;

Digits := 3

> vel;

v(t) = -37.5 + 37.5 e^(-.267 t)

As the rubber duck falls to the sea, time progresses.  Note what happens to v(t) in this process.

> for t from 2 by 6 to 38 do vel od;     t;    t := 't';    t;

v(t) = -15.5
v(8) = -36.6
v(14) = -36.6
v(20) = -37.4
v(26) = -37.5
v(32) = -37.5
v(38) = -37.5
44
t := t
t

Specifically,

> Limit( rhs(vel), t=infinity);        # try this without the above t := 't'

lim( -37.5 + 37.5 e^(-.267 t) , t = infinity)

> value(%);

-37.5

This is the terminal velocity of the rubber duck.   All objects falling to Earth only under the influence of gravity and a drag force either of the form D = cV or D = cV^2 will, given enough time, approach a terminal velocity.  Why?  Hint: Check out the following pictures.

duck force diagramduck force diagram
duck force diagramduck force diagram

Top

horzn_ln.gif (2407 bytes)

The earlier comment about when Newton's 2nd law creates run of the mill algebraic equations vs. when it creates differential equations is not entirely accurate; you can view all equations generated by Newton's 2nd law as differential equations.  For example, we'll use differential equations to find the position equation of the Kinematics 1 version of The Fall of the Rubber Duck, y(t) = 500 - (g t/2).

Without air resistance we have,

F[net] = - W

Or,

m diff( y(t), t, t ) = -m g

Or,

diff( v(t), t ) = -g

And this is the "simplest" form of a differential equation, which we used to derive the solution to the other types of differential equations.  Since it's an initial-value problem, it's solution is

Int( diff( v(s), s ), s=0..t ) = Int( -g, s=0..t )

v(t) - v(0) = -g t

v(t) = - g t

Now,

v(t) = diff( y(t), t )

Or,

diff( y(t), t ) = - g t

This is, of course, another instance of an initial-value problem in the simplest form.

Int( diff( y(s), s ), s=0..t ) = Int( -g, s=0..t )

y(t) - y(0) = - (g (t^2)/2)

Or,

y(t) = 500 - (g t^2)/2

Top

horzn_ln.gif (2407 bytes)

In the last set of commands we'll produce a multiple plot of the rubber duck's position vs. time graphs, both with and without air resistance.   It vividly illustrates the inaccuracy of problem solutions that fail to consider all contributing factors.

Note that, with air resistance,

> y(t)-y(0) = Int(rhs(subs(t=s,vel)),s=0..t);

y(t) - y(0) = Int( -37.5 + 37.5 e^(-.267 s), s=0..t )

> y_air := 500 + value(rhs(%));          # note that for t=0, y=500

y_air := 640 - 37.5 t - 140 e^(-.267 t)

And without air resistance,

> y_no_air := 500 - (g*t^2)/2;

y_no_air := 500 - 5 t^2

> y := plot({y_air,y_no_air}, t=0..17.5, labels=["t","y"], title="Height vs. Time"):

> with(plots):

> t1 := textplot([8.5,325,`with air resistance`], align={ABOVE,RIGHT}):

> t2 := textplot([12,-210,`without air resistance`], align={BELOW,LEFT}):

> display({y,t1,t2});

rubber duck's height vs. time multiplot

When considering the rubber duck's fall in a more realistic setting, the fall takes almost twice as long!  This says something about how often you can expect to see differential equations in any serious physics or engineering work that you do in the future.

 

horzn_ln.gif (2407 bytes)

Differential Equations Problem Set         Top         An Introduction to Diff Eqs

Homogeneous Linear Diff Eqs         Nonhomogeneous Linear Diff Eqs

Tutorial Index        Maple Index        How to...        Problem Set Index

Please send me any polite comments, suggestions, or corrections.