Chapter 4

Application of Mathematical Modeling

Introduction

How can we construct and use models in the mathematical world to help us better understand
real-world system? Before discussing how we link the two worlds together, let’s consider what
we mean by a real-world system and why we would be interested in constructing a mathematical
model for a system in the first place.

A system is an assemblage of objects joined in some regular interaction or independence. The

modeler is interested in understanding how a particular system works, what causes changes in
the system, and how sensitive the system is to certain changes.

In this section we will discuss different applications of modeling process to express real world
phenomenon.

4.1 Modeling using Proportionality

We introduced the concept of proportionality in chapter one in this section we use the concept of
proportionality in model formation.

Example1: Testing for Proportionality

Consider a spring-mass system, such as the one showing in figure 1.2. We conduct an
experiment to measure the stretch of the spring as a function of the mass (measured as weight)
placed on the spring. Consider the data collected for this experiment, dis playing in the table
below a scatter plot graph of the stretch or elongation of the spring versus the mass or weight
placed on it reveals an approximate straight line passing through the origin.

The data appear to follow the proportionality rule that elongation e is proportional to the mass m,
or symbolically, e  m. the straight line appears to pass through the origin.

1
Figure 4.1 spring-mass system

This geometric understanding allows us to look at the data to determine if proportionality is a

reasonable simplifying assumption and, if so, to estimate the slope k. In this case, the assumption
is valid, so we estimate the constant of proportionality by picking the two points (200, 3.25) and
(300, 4.875) as lying along the straight line. We calculate the slope of the line joining this points
4.875  3.25
as slope   0.01625
300  200

Figure 4.2Data from spring-mass system

Thus the constant of proportionality is approximately 0.0163 and we estimate our model as

e  0.0163 m

Now y  x if and only if y=kx for some constant k>0 (1)

2
Of course, if y  x , then x y because the constant k in equation (1) is greater than zero and then
1
x y . The following are other examples of proportionality relationships:
k

y  x 2 if and only if y  k1 x 2 for k1 a constant (2)

y  ln x if and only if y  k 2 ln x for k2 a constant (3)

y  e x if and only if y  k 3 e x for k3 a constant (4)

1
In equation (2), y  kx2 , k>0, so we also have x  y 1 2 because x  ( ) y 1 2 . This leads us to
k
consider how to link proportionalities together, a transitive rule for proportionality:

y  x and x  z, then y  z

Figure 4.3 Geometrical interpretation of y  x

Thus, any variables proportional to the same variables are proportional to one another.

Now let’s explore a geometric interpretation of proportionality. In equation (1), y  kx yields

k  y x . Thus, k may be interpreted as the tangent of the angle  depicted in Figure 4.3, and the
relation y  x defines a set of points along a line in the plane with angle of inclination  .

Comparing the general form of a proportionality relationship y  kx with the equation for a
straight line y  mx  b, we can see that the graph of a proportionality relationship is a line
(possibly extended) passing through the origin. If we plot the proportionality variables for
Models (2)-(4), we obtain the straight line graphs presented in Figure 4.4.

3
 Figure 4.4 Geometrical interpretations of models (a) (2), (b) (3), (c) (4)

Remark:-it is important to note that not just any straight line represents a proportionality
relationship: the y- intercept must be zero so that the line passes through the origin.

Example: Suppose we are interested in predicting the volume of water displaced by a boat as it
is loaded with cargo. Because a floating object displaces a volume of water equal to its weight,
we might be tempted to assume that the total volume y of displaced water is proportional to the
weight x of the added cargo. However, there is a flaw with that assumption because the unloaded
boat already displaces a volume of water equal to its weight. Although the graph of total volume
of displaced water versus weight of added cargo is given by a straight line, it is not given by a
line passing through the origin (Figure 4.5), so the proportionality assumption is incorrect.

Figure 4.5 A straight –line relationship between displaced volume and total weight, but it is not a
proportionality because the line fails to pass through the origin.

Example: Kepler’s Third Law

To assist in further understanding the idea of proportionality, let’s examine one of the famous
proportionalities from Table 1, Kepler’s third law. In 1601, the German astronomer Johannes
Kepler became director of the Prague Observatory.

4
Kepler’s had formulated his first two laws on the relative motion of the planet:

i) Each planet moves along an ellipse with the sun at one focus.

ii) For each planet, the line from the sun to the planet sweeps areas in equal times.

4.2 Modeling Using Geometric Similarity

Geometric similarity is a concept related to proportionality and can be useful to simplify the
mathematical modeling process.

Definition: - Two objects are said to be geometrically similar if there is a one to one
correspondence between points of the objects such that the ratio of distance between
corresponding points is constant for all possible pair’s points.

For example, consider the two boxes depicted in Figure 2.2.1. Let l denote the distance between
the points A and B in (a), and Let l’ be the distance between the corresponding points A’ and B’
in (b).Other corresponding points in the two figure, and associated distance between the points,
are marked the same way. For the boxes to be geometrically similar, it must be true that

l w h
   k. for some constant k>0
l ' w' h'

Figure4.6 Two geometrically similar objects

Let’s interpret the last result geometrically. In Figure 4.6, consider the triangles ABC and
A’B’C’. If the two boxes are geometrically similar, this triangle must be similar. The same
argument can be applied to any corresponding triangles, such as CBD and C’B’D’. Thus,
corresponding angles are equal for objects that are geometrically similar. In other words the

5
shape is the same for two geometrically similar objects, and one object is simply an enlarged
copy of the other.

One advantage that results when two objects are similar is a simplification in certain
computations, such as volume and surface area.

For example for the boxes in Figure 2.2.1, consider the following argument for the ratio of the
volumes V and V’:

V lwh
  k3 (5)
V ' l ' w' h'

Similarly, the ratio of their total surface areas S and S’ is given by

S 2lh  2wh  2wl

  k2 (6)
S ' 2l ' h'2w' h'2w' l '

Not only are these ratios immediately known once the scaling factor k has been specified but also
the surface area and volume may be expressed as proportionality in terms of some selected
characteristic dimensions. Let’s select the length l as the characteristic dimension. Then with
l l '  k, we have

S l2
 k2  2
S' l'

S S '2
Therefore,   cons tan t holds for any two geometrically similar objects. This is,
l 2 l '2
surface area is always proportional to the square of the characteristic dimension length:

S  l2

Likewise, volume is proportional to the length cubed:

V  l3

Thus, if we are interested in some function depending on an object’s length, surface area, and
volume, for example:

6
 y  f (l , S ,V )

We could express all the function arguments in terms of some selected characteristic
dimension, such as length, giving

y  (l , l 2 , l 3 )

Geometric similarity is a powerful simplifying assumption.

Example: Raindrops from a Motionless Cloud

Suppose we are interested in the terminal velocity of a raindrop from a motionless cloud.
Examining the free-body diagram, the only force acting on the raindrop are gravity and drag.

Assume that the atmospheric drag on the raindrop is proportional to its surface area S times the
square of its speed v. The mass m of the raindrop is proportional to the weight of the raindrop
(assuming constant gravity in Newton’s second law)

F  Fg  Fd  ma (7)

Under terminal velocity ( v  vt ), we have a=0 so equation (7) reduced to

Fg  Fd  0

Or

Fg  Fd

7
We are assuming that Fd  Sv2 and that Fg is proportional to weight w . Since m  w , we have

Fg  m.

Next we assume all the raindrops are geometrically similar. This assumption allows us to relate
area and volume so that

S  l 2 and V  l 3

For any characteristic dimension l. Thus, l  S 1 2  V 1 3 , which implies

S V 23

Because weight and mass are proportional to volume, the transitive rule for proportionality gives
S  m2 3

From the equation Fg  Fd , we now have m  m 2 3 vt . Solving for the terminal velocity, we have
2

1 1

m  v or m  vt
3 2
t
6

Therefore, the terminal velocity of the raindrop is proportional to its mass raised to the one-sixth
power.

4.3 Modeling Using Differential Equations

Introduction

Many phenomena can be described in a general way by saying that rates of change of the
endogenous variables depend on past and present values of the variables. These situations lead to
models involving differential and difference equations.

Having information relating a rate of change of a dependent variable with respect to one or more
independent variables and are interested in discovering the function relating the variables.

Example: If P represents the number of people in a large population with respect to time t then
it is reasonable to assume that the rate of change of the population with respect to time depends

8
on the current size of P as well as other factors like immigration, emigration, age, gender and
so on.

For ecological, economical, and other importance reasons, it is desirable to determine a relation
between P and t to make prediction about P.

If the present population size is denoted by P(t) and the population size at time t  t is
P (t  t ), then the change in population P during that time period t is given by

P  P (t  t )  P (t ) (1)

By assuming all other factors listed above neglected: P  P. and we can assume that during a
unit time period a certain percentage of the population reproduces while a certain percent age
dies. Suppose the constant of proportionality k is expressed as a percentage per unit time. Then
our proportionality assumption gives

P  P(t  t )  P(t )  kPt (2)

Equation (2) is a difference equation in which we are treating a discrete set of time period rather
than allowing t to vary continuously over some interval.

Assume that t does vary continuously so that we can take advantage of the calculus. Division of
equation (2) by t gives

P P(t  t )  P(t )
  kP (3)
t t

Next, allow t to approach zero. The definition of the derivatives gives the differential equation

P dP
lim   kP
t 0 t dt

where dP dt represents the instantaneous rate of change.

In modeling the derivatives is used in two distinct roles:

1. To represent the instantaneous rate of change in continuous problems

2. To approximate an average rate of change in discrete problems

9
The advantage of approximating an average rate of change by a derivative is that the calculus
often helps in uncovering a functional relationship between the variables under investigation.
The interpretation of the derivatives as an instantaneous rate of change is useful in many
modeling applications. The geometrical interpretation of the derivative as the slope of the line
tangent to the curve is useful for constructing numerical solutions.

4.3.1 The Derivative as a Rate of Change

The origin of the derivative lies in human kind’s curiosity about motion and our need to develop
a deeper understanding of motion. The search for the laws governing planetary motion, the study
of the pendulum and its application to clock building, and the law governing the flight of a
cannonball and so on.

Consider a particle whose distances from a fixed position depend on time t. Let the graph in
Figure 4.7 represent the distance s as a function of time t, and let (t1, s1) and (t2, s2) denotes two
points on the graph

Figure 4.7: Graph of distances as a function of time t

Define t  t 2  t1 , s  s2  s1 , and from the ratio s t Note that this ratio represents a rate:

an increment of distance traveled s over some increment of time t . That is, the ratio s t

represents the average velocity during the time period in question. Now the derivatives ds dt
evaluated at t=t1 is defined as

ds s
t t1  lim (4)
dt t 0 t

Using the interpretation of average velocity, we can see that at each state of using a smaller t
we are computing the average velocity over smaller and smaller intervals with left endpoints at t1
until, in the limit, we have the instantaneous velocity at t=t1. If we think of the motion of a

10
moving vehicle, this instantaneous velocity would correspond to the exact reading of its
speedometer at the instant t1.

1.3.2. The Derivatives as the Slope of the Tangent Line

Let’s consider another interpretation of the derivative. We consider s(t) simply as a curve. Let
examine a set of secant lines each emanating from the point A= (t1, s(t1)) on the curve. To each
secant there corresponds a pair of increments (t i , si ) as shown in Figure 4.8

Figure 4.8: The slope of each secant line approximates the slope of the tangent line to the curve
at the point A.

The lines AB, AC, and AD are secant lines. As t  0, these secant lines approach the line

tangent to the curve at the point A. Because the slope of each secant is s t , we may interpret
the derivative as the slope of the line tangent to the curve s(t) at the point A.

4.3.3: Some Mathematical Models Related to first order differential equations

I. Newton’s Law of Cooling

According to Newton’s empirical law of cooling, “the rate at which a body cools is proportional
to the difference between the temperature of the body and the temperature of the surrounding
medium.”

11
To model this physical law if we let T(t) represent the temperature of the body at any time t, T m
represent the constant temperature of the surrounding medium, and dT dt represent the rate at
which a body cools, then Newton’s law of cooling translates into the mathematical statement

dT dT
 T  Tm or  k (T  Tm ) (5)
dt dt

where k is a constant of proportionality. Since we have assumed the body is cooling, we must
have T  Tm , and so it stands to reason that k<0.

Example :- Cooling of a cake

When a cake is removed from an oven, its temperature is measured at 300 0 F . Three minute later
its temperature is 200 0 F . How long will it take for the cake to cool off to a room temperature of
70 0 F ?

Solution:- In (5) we make the identification Tm  70. We must then solve the initial-value
problem

dT
 k (T  70), T (0)  300 (6)
dt

And determine the value of k so that T (3)  200 .

Equation (6) is both linear and separable. Separating variables,

dT
 kdt
T  70

yields ln T  70  kt  c1 , and so T  70  c2 e kt . When t=0,T=300. So that 300  70  c2 gives

c2  230 and, therefore, T  70  230 e kt .

13 1 13
Finally, the measurement T (3)  200 . leads to e 3k  or k  ln  0.19018 .
23 3 23

Thus

12
 T (t )  70  230e 0.1918t . (7)

We note that (7) furnishes no finite solutions to T (t )  70. since lim T (t )  70.
t 

II. Population Growth

How many people will there be in the population of a certain country in n years? How many
births? In this section we build the simplest possible model for answering these questions.

Problem Identification: Suppose we know the population at some given time, for example, P0 at
time t=t0, and we are interested in predicting the population P at some future time t=t1. In other
words, we want to find a population function P(t) for t 0  t  t1 satisfying P(t 0 )  P0 .

Assumptions:

Consider some factors that pertain to population growth such as:

 the birthrate
 the death rate.

The birth rate and death rate are determined by different factors.

The birth rate is influenced by

 infant mortality rate,

 attitude toward and availability of contraceptives,
 attitudes toward abortion,
 health care during pregnancy, and so forth.

The death rate is affected by

 sanitations and public health,

 wars, pollutions,
 medicines,
 diet,
 psychological stress and anxiety, and so forth.

Other factors that influence population growth in a given region are:

13
  immigration and emigration,
 living space restrictions,
 availability of food and water, and
 epidemics.

For our model, let’s neglect all these latter factors. Now we will consider only the birthrate and
death rate. Because knowledge and technology have helped humankinds diminish the death rate
below the birthrate, human populations have tended to grow.

Let’s begin by assuming that a small unit time period a percentage b of the population is newly
born. And a percentage c of the population dies.

So the new population P (t  t ) is the old population P(t) plus the number of birth minus the
number of death during the time period t .

Symbolically

P(t  t )  P(t )  bP (t )t  cP (t )t

or

P
 bP  cP  (b  c) P  kP
t

From our assumptions the average rate of change of the population over an interval is
proportional to the size of the population.

Using the instantaneous rate of change to approximate the average rate of change, we have the
following differential equation model:

dP
 kP, P(t 0 )  P0 , t 0  t  t1 (8)
dt

where (for growth) k is a positive constant.

Solving the Model

14
We can separate the variables and rewrite equation (8) by moving all terms involving P and dP
to one side of the equation and all terms in t and dt to the other. This gives

dP
 kdt,
P

Integration of both sides of this last equation yields

ln P  kt  C (9)

For some constant C applying the condition P(t 0 )  P0 to equation (9) to find C results in

C  ln P0  kt0

Then, substitution for C in to equation (9) gives

ln P  kt  ln P0  kt

Or, simplifying algebraically,

P
ln  k (t  t 0 )
P0

Finally, we obtain the solution

P(t )  P0 e k (t t0 ) (10)

Equation (10), known as the Malthusian model of population growth, predicts that the
population grows exponentially with time.

Verifying the Model

P
Because ln  k (t  t 0 ) , our model predicts that if we plot ln P / P0 versus t  t 0 , a straight line
P0
passing through the origin with slope k should result. However, if we plot the population data for
the United States for several years, the model does not fit very well, especially in the later years.
In fact, the 1990 census for the population of the US was 248,710,000, and in 1970 it was
203,211,926. Substituting these values in equation (10) we can get the value of k as

15
 248,710,000
 e k (19901970)
203,211,926

Thus,

1 248,710,000
k ( ) ln  0.01
20 203,211,926

That is during the 20-year period from 1970-1990, population in the US was increasing at
average rate of 1%per year. We can use this information together with equation (10) to predict
the population for 2000, in this case t0=1990, P0=248,710,000 and k=0.01 yields

P(2000)  248,710,000e 0.01( 20001990)  303,775,080

The 2000 census for the population of the US was 281,400,000. Thus our prediction is off the
mark by approximately 8%. We can probably live with that magnitude of error, but let’s look
into the disaster future. Our model predicts that the population of the US will be 55,209 billion in
the year 2300, a population that far exceed current estimates of the maximum sustainable
population of the entire planet. Therefore, we are forced to conclude that our model is
unreasonable over the long term.

Refining the Model to Reflect Limited Growth

Let’s consider that the proportionality factor k, measuring the rate of population growth in
equation (8) is now no longer constant but a function of the population. As the population
increases and gets closer to the maximum population M, the rate k decreases. One simple sub-
model for k is the linear one

, where r is a constant.

Substitution in to equation (8) leads to

dP
 r (M  P) P, (11)*
dt

dP
Or  rdt (12)
P( M  P)

16
Again we assume the initial condition and Equation (11) is referred to logistic
growth.

By using partial fraction

1 1 1 1 
   
P( M  P) M  P M  P 

Thus, Equation (12) can be rewritten as

dP dP
  rMdt
P M P

which integrates to

ln P  ln M  P  rMt  C (13)

For some arbitrary constant C. Using the initial condition, we evaluate C in the case P<M:

P0
C  ln  rMt0
M P

Substituting in to Equation (13) and simplifying gives

P P0
ln  ln  rM (t  t 0 )
M P M P

P ( M  P0 )
Or ln  rM (t  t 0 )
P0 ( M  P )

Exponentiation both sides of this equation give

P ( M  P0 )
 e rM (t t0 )
P0 ( M  P )

Then,

P0 MerM (t t0 )  P( M  P0 )  P0 Pe rM (t t0 )

So that solving for the population P gives

17
 P0 MerM (t t0 )
P(t ) 
M  P0  P0 e rM (t t0 )

To estimate P as t  , we rewrite this last equation as

P0 M
P(t )  (14)
[ P0  ( M  P0 )e  rM ( t t0 ) ]

The graph of the limited growth Equation (14) is shown in Figure 4.9 for the case . Such a
curve is called a logistic curve.

Figure 4.9 Graph of the limited growth model.

III. Prescribing Drug Dosage

The problem of how much a drug dosage to prescribe and how often the dosage should be
administered is an important one in pharmacology. For most drugs there is a concentration below
which the drug is ineffective and a concentration above which the drug is dangerous.

Problem Identification:- How can the doses and the time between doses be adjusted to maintain
a safe but effective concentration of the drug in the blood.

The concentration in the blood resulting from a single dose of a drug normally decreases with
time as the drug is eliminated from the body (Figure 4.10).

18
Figure 4.10The concentration of a drug in the bloodstream decreases with time

We are interested in what happens to the concentration of the drug in the blood as doses are
given at regular intervals.

Now our aim is to Model this pharmacological idea in to Mathematical concept as follow

Let H denotes the highest save level of the drug in the bloodstream and

L denotes the lowest effective level of the drug in the bloodstream; it would be desirable to
prescribe a dosage C0 with time T between doses so that the concentration of the drug in the
blood stream remains between L and H over each dose period.

Let’s consider several ways in which the drugs might be administered. In Figure 4.11a the time
between doses is such that effectively there is no buildup of the drug in the system. On the other
hand in Figure 4.11b the interval between doses relative to the amount administered and the
decay rate of the concentration is such that a residual concentration exists at each time the drug is
taken.

19
Figure 4.11 Residual build up depends on the time interval between administrations of drug
doses

Our ultimate goal in prescribing drugs is to determine dose amounts and intervals between doses
and thereafter the concentration is maintained between L and H, as shown in Figure 4.12

Figure 4.12 Safe but effective levels of drug in the blood: C0 is the change in
concentration produced by one dose and T is the time interval between doses

Assumption:-To solve the problem we have identified, let’s consider the factors that determine
the concentration C(t) of the drug in the blood stream at any time t.

We begin with

C (t )  f (decayrate ,assimilati on rate , dosage amount , dosage int erval ,...)

And various other factors, including body weight and blood volume. To simplify our
assumptions, let’s assume body weight and blood volume are constants.

Next we determine sub models for decay rate and assimilation rate.

Sub-model for Decay Rate:-Consider the elimination of the drug from the bloodstream. This is
probably a discrete phenomenon, but let’s approximates it by a continuous function. Clinical

20
experiments have revealed that the decrease in the concentration of a drug in the blood stream
will be proportional to the concentration.

Mathematically this assumption means: C ' (t )  kC(t ) (15)

In this formula k is a positive constant called the elimination constant of the drug. Notice C ' (t )
is negative; it is to describe the decreasing concentration. In equation (15) the quantities
measured as follows: time t in hours, C(t) is (mg/ml), C ' (t ) is mg/[Link], and k is hr-1.

Assume that the concentration H and L can be determined experimentally for a given population,
such as an age group. Then set the drug concentration for a single dose at the level

C0  H  L (16)

If we assume that C0 is the concentration at t=0, then we have the model

dC
 kC, C (0)  C0 (17)
dt

The variables can be separable in Equation (17). The solution of the model gives

C (t )  C0 e  kt (18)

The graph of C(t) looks like the one in Figure 4.13

Figure 4.13 Exponential model for decay of drug concentration with time

Sub-model for Assimilation Rate:-Having made an assumption about how drug concentration
decreases with time, let’s consider how they increase again when drugs are administered. Our
initial assumption is that when a drug is taken, it is diffused rapidly throughout the blood that the
graph looks vertical. That is, we assume an instantaneous rise in concentration whenever a drug
is administered.

21
Now let’s see how the drug accumulates in the bloodstream with repeated doses.

Drug Accumulation with Repeated Doses:- Consider what happens to the concentration C(t)
when a dose that is capable of raising the concentration by C0 mg/ml each time it is given is
administered regularly at fixed time intervals of length T.

Suppose at time t=0 the first dose is administered. According to model (18),

After T hours have elapsed, the residual R1  C0 e  kT remains in the blood, and then the
second dose is administered. Because of our assumption concerning the increase in drug
concentration, the level jumps to C1  C0  C0 e  kT then after T hours elapse again, the residual

R2  C1e  kT  C0 e  kT  C0 e 2kT remains in the blood. In similar fashion this continues for n

successive time we determine a formula for the nth residual Rn.

Rn  C 0 e  kT  C 0 e 2 kT  ...  C 0 e  nkT
(19)
 C 0 e kT (1  r  r 2  ...  r n 1 )

Where r  e  kT .

Figure 4.14 One possible effect of repeating equal doses.

Algebraically,

n 1 1 rn
1  r  r  ...  r
2

1 r

So substitution for r in Equation (19)gives the result

C 0 e  kT (1  e  nkT )
Rn  (20)
1  e kT

22
When n   , the number e  nkT is close to 0. As a result, the sequence of Rn’s has a limiting
value, which we call R:

C0 e  kT
R  lim Rn 
n  1  e kT

C0
or R (21)
e 1
kT

In summary, if a dose is capable of raising the concentration by C0 mg/ml is repeated at intervals

of T hours, then the limiting value R of the residual concentration is given by Equation (21)

Determining the Dose Schedule:- The concentration Cn-1 at the beginning of the nth interval is
given by

Cn1  C0  Rn1 (22)

Cn-1 to approach H as n becomes large. That is,

H  lim Cn1  lim (C0  Rn1 )  C0  R

n n

Combining this last result with C0  H  L yields

RL (23)

A meaningful way to examine what happens to the residual concentration R for different
intervals T between doses is to examine R in comparison with C0 , the change in concentration to
each dose. To make this comparison, we form the dimensionless ratio

R 1
 kT (24)
C0 e  1

Then substitution of R  L and C0  H  L in Equation (21) yields

23
 H L
L
e kT  1

We then solve the preceding equation for e kT to obtain

e kT  H / L

Taking the logarithm of both sides of this last equation and dividing the result by k gives the
desired dose schedule:

1 H
T ln (25)
k L

To reach an effective level rapidly, administer a dose, often called a loading dose, which will
immediately produce a blood concentration of H mg/ml. This medication can be followed every
1 H
T ln hour by a dose that raises the concentration by C0  H  L mg/ml.
k L

4.4 Modeling with Higher-order differential Equations

In this section we are going to consider several linear dynamical systems in which each
mathematical model is a second order differential equation with constant coefficients

d2y dy
a2 2
 a1  a0 y  g (t ).
dt dt

Recall that the function g is the input or forcing function of the system. A solution of the
differential equation on an interval containing t0 and satisfying prescribed initial conditions
y(t0)=y0, y’(t0)=y1 is the output or response of the system.

4.4.1 Spring/Mass Systems: Free Un-damped Motion

Hooke's Law: Suppose, as in Figure 4.15, that a mass m1 is attached to a flexible spring
suspended from a rigid support. When m1 is replaced with a different mass m2, the amount of
stretch, or elongation, of the spring will of course be different.

24
 Figure 4.15: Spring/Mass Systems

By Hooke's law the spring itself exerts a restoring force F opposite to the direction of elongation
and proportional to the amount s of elongation.

Simply stated, F  ks, where k is a constant of proportionality called the spring constant.
Although masses with different weights stretch a spring by different amounts, the spring is
essentially characterized by the number k. For example, if a mass weighing 10kg stretches a
spring 1 / 2 cm, then 10 = k(1/2) implies k = 20 kg/cm. Necessarily then, a mass weighing, say, 8
kg stretches the same spring 2/5cm.

Newton's Second Law After a mass m is attached to a spring; it stretches the spring by an
amount s and attains a position of equilibrium at which its weight W is balanced by the restoring
force ks. Recall that weight is defined by W = mg, where mass is measured in slugs, kilograms,or
grams and g = 32 ft/s 2 , 9.8 m/s 2 , or 980 cm/s 2 , respectively. As indicated in Figure 2.4.2(b),
the condition of equilibrium is mg  ks. or mg  ks  0. If the mass is displaced by an amount x
from its equilibrium position, the restoring force of the spring is then k(x + s). Assuming that
there are no retarding forces acting on the system and assuming that the mass vibrates free of
other external forces-free motion-we can equate Newton's second law with the net, or resultant,
force of the restoring force and the weight:

d 2x
m  k ( s  x)  mg  kx  mg  ks  kx (1)
dt 2 
 
zero

25
The negative sign in (1) indicates that the restoring force of the spring acts opposite to the
direction of motion. Furthermore, we can adopt the convention that displacements measured
below the equilibrium position are positive. See Figure 4.16

Figure 4.16 Mass spring system

Differential Equation of Free Undamped Motion: By dividing (1) by the mass m we obtain
the second-order differential equation d 2 x dt 2  (k / m) x  0 or

d 2x
2
2x  0 (2)
dt

Where  2  k / m. . Equation (2) is said to describe simple harmonic motion or free undamped
motion. Two obvious initial conditions associated with (2) are x (0)   , the amount of initial
displacement, and x' (0)   , the initial velocity of the mass. For example, if   0,   0, the
mass starts from a point below the equilibrium position with an imparted upward velocity. If
  0,   0, the mass is released from rest from a point  units above the equilibrium

position, and so on.

Solution and Equation of Motion: To solve equation (2) we note that the solutions of the
auxiliary equation m 2   2  0 are the complex numbers m1  i, m2  i . We find the general
solution of (2) to be

26
 x(t )  c1 cost  c2 sin t. (3)

The period of free vibrations described by (3) is T  2 /  , and the frequency is

f  1 / T   / 2 . For example, for x(t )  2 cos 3t  4 sin 3t. the period is 2 / 3 and the

frequency is 3 / 2 . . The former number means that the graph of x(t) repeats every 2 / 3 units;
the latter number means that there are 3 cycles of the graph every 2 units or, equivalently, that

the mass undergoes 3 / 2 . complete vibrations per unit time. In addition, it can be shown that the
period 2 / 3 is the time interval between two successive maxima of x(t). Keep in mind that a
maximum of x(t) is a positive displacement corresponding to the mass's attaining a maximum
distance below the equilibrium position, whereas a minimum of x(t) is a negative displacement
corresponding to the mass's attaining a maximum height above the equilibrium position. We
refer to either case as an extreme displacement of the mass. Finally, when the initial conditions
are used to determine the constants c1 and c2 in (3), we say that the resulting particular solution
or response is the equation of motion.

Example 1 Interpretation of an IVP

Solve and interpret the initial-value problem

d 2x
 16 x  0, x(0)  10, x' (0)  0.
dt 2

SOLUTION The problem is equivalent to pulling a mass on a spring down 10 units below the
equilibrium position, holding it until t=0, and then releasing it from rest. Applying the initial
conditions to the solution

x(t )  c1 cos 4t  c2 sin 4t.

gives x(0)  c1.1  c2 .0 so that c1=10. Hence

x(t )  10 cos 4t  c2 sin 4t

From x' (t )  40 sin 4t  4c2 cos 4t we see that x' (0)  0  4c2 .1 and so c2=0. Therefore the
equation of motion is x(t )  10 cos 4t.

27
The solution clearly shows that once the system is set in motion, it stays in motion, with the mass
bouncing back and forth 10 units on either side of the equilibrium position x=0. As shown in
Figure 4.17 the period of oscillation is 2 / 4   / 2 s.

Figure 4.17

EXAMPLE2 Free Un-damped Motion

A mass weighing 2 pounds stretches a spring 6 inches. At t= 0 the mass is released from a point
8 inches below the equilibrium position with an upward velocity of 4/3ft/s. Determine the
equation of free motion.

SOLUTION Because we are using the engineering system of units, the measurements given in
terms of inches must be converted into feet: 6 in. = 1/2ft; 8 in. = 2/3ft. In addition, we must
convert the units of weight given in pounds into units of mass. From m = W/g we have
2 1
m  slug. Also, from Hooke's law, 2 = k(1/2) implies that the spring constant is k = 4
32 16
lb/ft. Hence (1) gives

1 d 2x d 2x
 4 x or  64 x  0
16 dt 2 dt 2

2 4
The initial displacement and initial velocity are x(0)  , x' (0)   , where the negative sign in
3 3
the last condition is a consequence of the act that the mass is given an initial velocity in the
negative, or upward, direction.

Now  2  64 or   16 , so that the general solution of the differential equation is

28
 x(t )  c1 cos8t  c2 sin 8t. (4)

2 1
Applying the initial conditions to x(t) and x' (t) gives c1  and c2   . Thus the equation of
3 6
motion is

2 1
x(t )  cos 8t  sin 8t. (5)
3 6
Alternative Form of x(t) When, c1  0 and c2  0, the actual amplitude A of free vibrations is
not obvious from inspection of Equation(3). For example, although the mass in Example 2 is
initially displaced 2/3 foot beyond the equilibrium position, the amplitude of vibrations is a
number larger than 2/3. Hence it is often convenient to convert a solution of form (3) to the
simpler form

x(t )  A sin( t   ), (6)

where A  c1  c2 and  is a phase angle defined by

2 2

c1 
sin  
A  c1
 tan   , (7)
c c2
cos   2 
A

To verify this we expand (6) by the addition formula for the sine function:

A sin t cos   A cos t sin   ( A sin  ) cos t  ( A cos  ) sin t. (8)

It follows from Figure 2.4.5 that if  is defined by

c1 c1 c2 c2
sin    , cos    , then (8) becomes
c1  c 2 A c1  c 2 A
2 2 2 2

c1 c
A cos t  A 2 sin t  c1 cos t  c2 sin t  x(t ).
A A

29
4.4.2 Spring/Mass Systems: Free Damped Motion
The concept of free harmonic motion is somewhat unrealistic since the motion described by
equation (1) assumes that there are no retarding forces acting on the moving mass. Unless the
mass is suspended in a perfect vacuum, there will be at least a resisting force due to the
surrounding medium. As Figure 4.18 shows, the mass could be suspended in a viscous medium
or connected to a dashpot damping device.

Differential Equation of Free Damped Motion In the study of mechanics, damping forces
acting on a body are considered to be proportional to a power of the instantaneous velocity. In
particular, we shall assume throughout the subsequent discussion that this force is given by a
constant multiple of dx / dt. When no other external forces are impressed on the system, it follows
from Newton's second law that

d 2x dx
m 2  kx   , (9)
dt dt

where  is a positive damping constant and the negative sign is a consequence of the fact that
the damping force acts in a direction opposite to the motion.

Dividing (10) by the mass m, we find the differential equation of free damped motion is
d 2 x / dt 2  ( / m)dx / dt  (k / m) x  0. or

d 2x dx
2
 2   2 x  0, (10)
dt dt

 k
Where 2  , 2  (11)
m m

The symbol 2 is used only for algebraic convenience since the auxiliary equation is

m 2  2m   2  0 and the corresponding roots are then m1    2   2 ,

m 2     2   2 , .

We can now distinguish three possible cases depending on the algebraic sign of 2   2 . Since
each solution contains the damping factor e t ,   0. the displacements of the mass become
negligible for large time.

30
CASE I: 2   2  0 In this situation the system is said to be over-damped since the damping
coefficient  is large when compared to the spring constant k. The corresponding solution of

(10) is x(t )  c1e m1t  c2 e m2t or

x(t )  e t (c1e 2  2t

 c2 e  2  2t
). (12)

This equation represents a smooth and non oscillatory motion. Figure 4.19 shows two possible
graphs of x(t).

Figure 4.18 mass suspended in a viscous medium

Figure 4.19 two possible graphs of x(t) in Equation (12)

Case II: 2   2  0 The system is said to be critically damped since any slight decrease in the
damping force would result in oscillatory motion. The general solution of (10) is
x(t )  c1e m1t  c2te m1t or

x(t )  e t (c1  c2 t ). (13)

31
Notice that the motion is quite similar to that of an over damped system. It is also apparent from
(13) that the mass can pass through the equilibrium position at most one time.

CASE III: 2   2  0 In this case the system is said to be under- damped since the damping
coefficient is small compared to the spring constant. The roots ml and m2 are now complex:

m1     2  2 i , m1     2  2 i .

Thus the general solution of Equation (10) is

x(t )  e  t (c1 cos  2  2 t  c 2 sin  2  2 t ). (14)

As indicated in Figure 2.4.8 the motion described by (15) is oscillatory; but because of the
coefficient e  t , the amplitudes of vibration  0 as t  .

Figure 4.20

EXAMPLE 4 Over-damped Motions

It is readily verified that the solution of the initial-value problem

d 2x dx
2
 5  4 x  0, x(0)  1, x' (0)  1
dt dt

5 2
x(t )  e t  e 4t (15)
3 3

The problem can be interpreted as representing the over-damped motion of a mass on a spring.
The mass starts from a position 1 unit below the equilibrium position with a downward velocity
of 1 ft/s.

32
To graph x(t) we find the value of t for which the function has an extremum-that is, the value of
time for which the first derivative (velocity) is zero. Differentiating (15) gives
5 8 8 1 8
x' (t )   e t  e 4t so that x' (t )  0 implies e 3t  or t  ln  0.157 . It follows from the
3 3 5 3 5
first derivative test, as well as our physical intuition, that x(0.157)=1.069 ft is actually a
maximum. In other words, the mass attains an extreme displacement of 1.069 feet below the
equilibrium position.

EXAMPLE 5 Critically Damped Motions

An 8-pound weight stretches a spring 2 feet. Assuming that a damping force numerically equal to
2 times the instantaneous velocity acts on the system, model the equation of motion if the weight
is released from the equilibrium position with an upward velocity of 3 ft/s.

SOLUTION From Hooke's law we see that 8=k(2)gives k=4 lb/ft and that W=mg gives
8 1
m  slugs . The differential equation of motion is then
32 4

1 d 2x dx d 2x dx
2
 4 x  2 or 2  8  16 x  0 (16)
4 dt dt dt dt

The auxiliary equation for (16) is m 2  8m  16  (m  4) 2  0 so that m1  m2  4 . Hence the

system is critically damped and

x(t )  c1e 4t  c2 te 4t . (17)

Applying the initial conditions x(0)=0 and x’(0)=-3, we find, in turn, that Cl =0 and C2=- 3. Thus
the equation of motion is modeled as

x(t )  3te 4t . (18)

Example 6 Under-damped Motion

A 16-pound weight is attached to a 5-foot-Iong spring. At equilibrium the spring measures 8.2
feet. If the weight is pushed up and released from rest at a point 2 feet above the equilibrium

33
position, find the displacements x(t) if it is further known that the surrounding medium offers a
resistance numerically equal to the instantaneous velocity.

SOLUTION The elongation of the spring after the weight is attached is 8.2-5=3.2ft, so it
follows from Hooke's law that 16=k(3.2) or k =5 lb/ft. In addition, m=16/32=1/2 slug so that the
differential equation is given by

1 d 2x dx d 2x dx
2
 5 x  or 2
 2  10 x  0 (19)
2 dt dt dt dt

Proceeding, we find that the roots of m 2  2m  16  0 are m1  1  3i and m2  1  3i

which then implies the system is under damped and

x(t )  e t (c1 cos 3t  c2 sin 3t ). (20)

Finally, the initial conditions x(0)=-2 and x’(0)=0 yields c1=-2 and c2=-2/3, so the equation of
motion is

2
x(t )  e t (2 cos 3t  sin 3t ) (21)
3

4.4.3 Spring/Mass Systems: Driven Motion

Differential Equation of Driven Motion with Damping Suppose we now take into
consideration an external force f(t) acting on a vibrating mass on a spring. For example, f(t)
could represent a driving force causing an oscillatory vertical motion of the support of the spring.
See Figure 4.18. The inclusion of f(t) in the formulation of Newton's second law gives the
differential equation of driven or forced motion:

d 2x dx
m 2
 kx    f (t ) (22)
dt dt

Dividing (24) by m gives

d 2x dx
2
 2   2 x  F (t ) (23)
dt dt

where F(t) == f(t)/m and, as in the preceding section, 2   / m,  2  k / m .

34
To solve the latter non homogeneous equation we can use either the method of undetermined
coefficients or variation of parameters.

EXAMPLE 7 Interpretation of an Initial-Value Problem

Interpret and solve the initial-value problem I

1 d 2x dx 1
2
 1.2  2 x  5 cos 4t , x(0)  , x' (0)  1 (24)
5 dt dt 2

SOLUTION We can interpret the problem to represent a vibrational system consisting of a mass
(m=1/5slug or kilogram) attached to a spring (k = 2 lb/ft or N/m). The mass is released from
rest1/2 unit (foot or meter) below the equilibrium position. The motion is damped (   1.2) and
is being driven by. An external periodic ( T   / 2 s ) force beginning at t = 0. Intuitively we
would expect that even with damping the system would remain in motion until such time as the
forcing function was "turned off," in which case the amplitudes would diminish.

However, as the problem is given, f (t )  5 cos 4t will remain "on" forever.

We first multiply the differential equation in (24) by 5 and solve

d 2x dx
2
 6  10 x  0,
dt dt

by the usual methods. Since m1  3  i and m2  3  i , it follows that

xc (t )  e 3t (c1 cos t  c2 sin t ).

Using the method of undetermined coefficients, we assume a particular solution of the form
x p (t )  A cos 4t  B sin 4t ). Now

x' p (t )  4 A sin 4t  4B cos 4t , x' ' p (t )  16 A cos 4t  16 B sin 4t.

so that x p ' '6 x p '10 x p  (6 A  24 B) cos 4t  (24 A  6B) sin 4t  25 cos 4t

The resulting system of equations -6A + 24B = 25, - 24A - 6B = 0

yields A = -25/102 and B = 50/51 . It follows that

35
 25 50
x (t )  e 3t (c1 cos t  c2 sin t )  cos 4t  sin 4t
102 51

38
When we set t = 0 in the above equation, we obtain c1  . By differentiating the expression
51
86
and then setting t = 0, we also find that c2  
51

38 86 25 50
Therefore the equation of motion is x (t )  e 3t ( cos t  sin t )  cos 4t  sin 4t
51 51 102 51

36