Systems of Differential Equations
11.1: Examples of Systems
11.2: Basic First-order System Methods 11.3: Structure of Linear Systems 11.4: Matrix Exponential
11.5: The Eigenanalysis Method for x′ =Ax 11.6: Jordan Form and Eigenanalysis
11.7: Nonhomogeneous Linear Systems 11.8: Second-order Systems
11.9: Numerical Methods for Systems
Linear systems. A linear system is a system of differential equa- tions of the form
x′1 = a11x1 + · · · + a1nxn + f1, x′2 = a21x1 + · · · + a2nxn + f2,
... ... · · · ... ... x′m = am1x1 + · · · + amnxn + fm, (1)
where′ =d/dt. Given are the functionsaij(t) andfj(t) on some interval a < t < b. The unknowns are the functionsx1(t), . . . , xn(t).
The system is called homogeneous if all fj = 0, otherwise it is called non-homogeneous.
Matrix Notation for Systems. A non-homogeneous system of linear equations (1) is written as the equivalent vector-matrix system
x′ =A(t)x+f(t), where
x=
x1
... xn
, f =
f1
... fn
, A=
a11 · · · a1n ... · · · ... am1 · · · amn
.
11.1 Examples of Systems
Brine Tank Cascade ... 521
Cascades and Compartment Analysis ... 522
Recycled Brine Tank Cascade ... 523
Pond Pollution ... 524
Home Heating ... 526
Chemostats and Microorganism Culturing ... 528
Irregular Heartbeats and Lidocaine ... 529
Nutrient Flow in an Aquarium ... 530
Biomass Transfer ... 531
Pesticides in Soil and Trees ... 532
Forecasting Prices ... 533
Coupled Spring-Mass Systems ... 534
Boxcars ... 535
Electrical Network I ... 536
Electrical Network II ... 537
Logging Timber by Helicopter ... 538
Earthquake Effects on Buildings ... 539
Brine Tank Cascade
Let brine tanks A,B,C be given of volumes 20, 40, 60, respectively, as in Figure 1.
water
C A
B
Figure 1. Three brine tanks in cascade.
It is supposed that fluid enters tank A at rate r, drains from A to B at rate r, drains from B to C at rate r, then drains from tank C at rate r. Hence the volumes of the tanks remain constant. Let r = 10, to illustrate the ideas.
Uniform stirring of each tank is assumed, which implies uniform salt concentration throughout each tank.
Letx1(t),x2(t),x3(t) denote the amount of salt at timet in each tank.
We suppose added to tank A water containing no salt. Therefore, the salt in all the tanks is eventually lost from the drains. The cascade is modeled by thechemical balance law
rate of change = input rate − output rate.
Application of the balance law, justified below incompartment analysis, results in the triangular differential system
x′1 =−1 2x1, x′2 =1
2x1− 1 4x2, x′3= 1
4x2− 1 6x3.
The solution, to be justified later in this chapter, is given by the equations x1(t) =x1(0)e−t/2,
x2(t) =−2x1(0)e−t/2+ (x2(0) + 2x1(0))e−t/4, x3(t) = 3
2x1(0)e−t/2−3(x2(0) + 2x1(0))e−t/4 + (x3(0)−3
2x1(0) + 3(x2(0) + 2x1(0)))e−t/6.
Cascades and Compartment Analysis
A linear cascadeis a diagram of compartments in which input and output rates have been assigned from one or more different compart- ments. The diagram is a succinct way to summarize and document the various rates.
The method of compartment analysis translates the diagram into a system of linear differential equations. The method has been used to derive applied models in diverse topics like ecology, chemistry, heating and cooling, kinetics, mechanics and electricity.
The method. Refer to Figure 2. A compartment diagram consists of the following components.
Variable Names Eachcompartment is labelled with a variableX.
Arrows Each arrow is labelled with aflow rate R.
Input Rate An arrowhead pointing at compartment X docu- mentsinput rate R.
Output Rate An arrowhead pointing away from compartmentX documentsoutput rateR.
0
x3
x2 x1
x3/6 x2/4 x1/2
Figure 2. Compartment analysis diagram.
The diagram represents the classical brine tank problem of Figure 1.
Assembly of the single linear differential equation for a diagram com- partment X is done by writingdX/dt for the left side of the differential equation and then algebraically adding the input and output rates to ob- tain the right side of the differential equation, according to the balance
law dX
dt = sum of input rates−sum of output rates
By convention, a compartment with no arriving arrowhead has input zero, and a compartment with no exiting arrowhead has output zero.
Applying the balance law to Figure 2 gives one differential equation for each of the three compartments x1 , x2 , x3 .
x′1= 0− 1 2x1, x′2= 1
2x1− 1 4x2, x′3= 1
4x2− 1 6x3.
Recycled Brine Tank Cascade
Let brine tanks A,B,C be given of volumes 60, 30, 60, respectively, as in Figure 3.
A
B
C
Figure 3. Three brine tanks in cascade with recycling.
Suppose that fluid drains from tank A to B at rate r, drains from tank B to C at rate r, then drains from tank C to A at rate r. The tank volumes remain constant due to constant recycling of fluid. For purposes of illustration, let r= 10.
Uniform stirring of each tank is assumed, which implies uniform salt concentration throughout each tank.
Let x1(t),x2(t),x3(t) denote the amount of salt at time tin each tank.
No salt is lost from the system, due to recycling. Using compartment
analysis, the recycled cascade is modeled by the non-triangular system x′1 = −1
6x1 + 1
6x3, x′2 = 1
6x1 − 1 3x2,
x′3 = 1
3x2 − 1 6x3. The solution is given by the equations
x1(t) =c1+ (c2−2c3)e−t/3cos(t/6) + (2c2+c3)e−t/3sin(t/6), x2(t) = 1
2c1+ (−2c2−c3)e−t/3cos(t/6) + (c2−2c3)e−t/3sin(t/6), x3(t) =c1+ (c2+ 3c3)e−t/3cos(t/6) + (−3c2+c3)e−t/3sin(t/6).
At infinity, x1 = x3 = c1, x2 = c1/2. The meaning is that the total amount of salt is uniformly distributed in the tanks, in the ratio 2 : 1 : 2.
Pond Pollution
Consider three ponds connected by streams, as in Figure 4. The first pond has a pollution source, which spreads via the connecting streams to the other ponds. The plan is to determine the amount of pollutant in each pond.
1 3 2
f(t)
Figure 4. Three ponds 1, 2, 3 of volumes V1,V2,V3 connected by streams. The pollution source f(t)is in pond 1.
Assume the following.
• Symbol f(t) is the pollutant flow rate into pond 1 (lb/min).
• Symbols f1,f2,f3 denote the pollutant flow rates out of ponds 1, 2, 3, respectively (gal/min). It is assumed that the pollutant is well-mixed in each pond.
• The three ponds have volumes V1,V2,V3 (gal), which remain con- stant.
• Symbolsx1(t),x2(t),x3(t) denote the amount (lbs) of pollutant in ponds 1, 2, 3, respectively.
The pollutant flux is the flow rate times the pollutant concentration, e.g., pond 1 is emptied with fluxf1 timesx1(t)/V1. A compartment analysis is summarized in the following diagram.
x2
x3
x1 f1x1/V1 f(t)
f3x3/V3 f2x2/V2
Figure 5. Pond diagram.
The compartment diagram represents the three-pond pollution problem of Figure 4.
The diagram plus compartment analysis gives the following differential equations.
x′1(t) = f3
V3 x3(t)− f1
V1x1(t) +f(t), x′2(t) = f1
V1
x1(t)− f2 V2
x2(t), x′3(t) = f2
V2 x2(t)− f3 V3x3(t).
For a specific numerical example, take fi/Vi = 0.001, 1 ≤ i ≤ 3, and let f(t) = 0.125 lb/min for the first 48 hours (2880 minutes), thereafter f(t) = 0. We expect due to uniform mixing that after a long time there will be (0.125)(2880) = 360 pounds of pollutant uniformly deposited, which is 120 pounds per pond.
Initially, x1(0) = x2(0) = x3(0) = 0, if the ponds were pristine. The specialized problem for the first 48 hours is
x′1(t) = 0.001x3(t)−0.001x1(t) + 0.125, x′2(t) = 0.001x1(t)−0.001x2(t),
x′3(t) = 0.001x2(t)−0.001x3(t), x1(0) =x2(0) =x3(0) = 0.
The solution to this system is x1(t) =e−20003t 125√
3 9 sin
√3t 2000
!
−125 3 cos
√3t 2000
!!
+125 3 + t
24, x2(t) =−250√
3
9 e−20003t sin(
√3t 2000) + t
24, x3(t) =e−20003t 125
3 cos
√3t 2000
!
+125√ 3 9 sin
√3t 2000
!!
+ t
24 −125 3 . After 48 hours elapse, the approximate pollutant amounts in pounds are
x1(2880) = 162.30, x2(2880) = 119.61, x3(2880) = 78.08.
It should be remarked that the system above is altered by replacing 0.125 by zero, in order to predict the state of the ponds after 48 hours. The
corresponding homogeneous system has an equilibrium solutionx1(t) = x2(t) = x3(t) = 120. This constant solution is the limit at infinity of the solution to the homogeneous system, using the initial valuesx1(0)≈ 162.30, x2(0)≈119.61, x3(0)≈78.08.
Home Heating
Consider a typical home with attic, basement and insulated main floor.
Attic MainFloor Basement
Figure 6. Typical home with attic and basement.
The below-grade basement and the attic are un-insulated.
Only the main living area is insulated.
It is usual to surround the main living area with insulation, but the attic area has walls and ceiling without insulation. The walls and floor in the basement are insulated by earth. The basement ceiling is insulated by air space in the joists, a layer of flooring on the main floor and a layer of drywall in the basement. We will analyze the changing temperatures in the three levels using Newton’s cooling law and the variables
z(t) = Temperature in the attic,
y(t) = Temperature in the main living area, x(t) = Temperature in the basement, t= Time in hours.
Initial data. Assume it is winter time and the outside temperature in constantly 35◦F during the day. Also assumed is a basement earth temperature of 45◦F. Initially, the heat is off for several days. The initial values at noon (t= 0) are thenx(0) = 45, y(0) =z(0) = 35.
Portable heater. A small electric heater is turned on at noon, with thermostat set for 100◦F. When the heater is running, it provides a 20◦F rise per hour, therefore it takes some time to reach 100◦F (probably never!). Newton’s cooling law
Temperature rate = k(Temperature difference)
will be applied to five boundary surfaces: (0) the basement walls and floor, (1) the basement ceiling, (2) the main floor walls, (3) the main
floor ceiling, and (4) the attic walls and ceiling. Newton’s cooling law gives positive cooling constants k0,k1,k2,k3,k4 and the equations
x′ = k0(45−x) +k1(y−x),
y′ = k1(x−y) +k2(35−y) +k3(z−y) + 20, z′ = k3(y−z) +k4(35−z).
The insulation constants will be defined ask0 = 1/2,k1 = 1/2,k2 = 1/4, k3 = 1/4, k4 = 1/2 to reflect insulation quality. The reciprocal 1/k is approximately the amount of time in hours required for 63% of the temperature difference to be exchanged. For instance, 4 hours elapse for the main floor. The model:
x′ = 1
2(45−x) +1
2(y−x), y′ = 1
2(x−y) +1
4(35−y) +1
4(z−y) + 20, z′ = 1
4(y−z) +1
2(35−z).
The homogeneous solution in vector form is given in terms of constants a = (7−√
21)/8, b = (7 +√
21)/8, c = −1 +√
21, d =−1−√ 21 and arbitrary constants c1,c2,c3 by the formula
xh(t) yh(t) zh(t)
=c1e−t/2
1 2 2
+c2e−at
4 c 2
+c3e−bt
4 d 2
. A particular solution is an equilibrium solution
xp(t) yp(t) zp(t)
=
455 2758 1854 4
.
The homogeneous solution has limit zero at infinity, hence the temper- atures of the three spaces hover around x= 57, y= 69, z= 46 degrees Fahrenheit. Specific information can be gathered by solving forc1,c2,c3
according to the initial data x(0) = 45, y(0) =z(0) = 35. The answers are
c1=−85
24, c2 =−25
24− 115√ 21
168 , c3=−25
24 +115√ 21 168 . Underpowered heater. To the main floor each hour is added 20◦F, but the heat escapes at a substantial rate, so that after one houry ≈51◦F.
After five hours, y ≈65◦F. The heater in this example is so inadequate that even after many hours, the main living area is still under 69◦F.
Forced air furnace. Replacing the space heater by a normal furnace adds the difficulty of switches in the input, namely, the thermostat
turns off the furnace when the main floor temperature reaches 70◦F, and it turns it on again after a 4◦F temperature drop. We will suppose that the furnace has four times the BTU rating of the space heater, which translates to an 80◦F temperature rise per hour. The study of the forced air furnace requires two differential equations, one with 20 replaced by 80 (DE 1, furnace on) and the other with 20 replaced by 0 (DE 2, furnace off). The plan is to use the first differential equation on time interval 0≤t≤t1, then switch to the second differential equation for time interval t1 ≤ t ≤ t2. The time intervals are selected so that y(t1) = 70 (the thermostat setting) and y(t2) = 66 (thermostat setting less 4 degrees). Numerical work gives the following results.
Time in minutes Main floor temperature Model Furnace
31.6 70 DE 1 on
40.9 66 DE 2 off
45.3 70 DE 1 on
55.0 66 DE 2 off
The reason for the non-uniform times between furnace cycles can be seen from the model. Each time the furnace cycles, heat enters the main floor, then escapes through the other two levels. Consequently, the initial conditions applied to models 1 and 2 are changing, resulting in different solutions to the models on each switch.
Chemostats and Microorganism Culturing
A vessel into which nutrients are pumped, to feed a microorganism, is called a chemostat1. Uniform distributions of microorganisms and nutrients are assumed, for example, due to stirring effects. The pumping is matched by draining to keep the volume constant.
In a typical chemostat, one nutrient is kept in short supply while all others are abundant. We consider here the question of survival of the organism subject to the limited resource. The problem is quantified as follows:
x(t) =the concentration of the limited nutrient in the vessel, y(t) =the concentration of organisms in the vessel.
1The October 14, 2004 issue of the journal Nature featured a study of the co- evolution of a common type of bacteria, Escherichia coli, and a virus that infects it, called bacteriophage T7. Postdoctoral researcher Samantha Forde set up ”micro- bial communities of bacteria and viruses with different nutrient levels in a series of chemostats – glass culture tubes that provide nutrients and oxygen and siphon off wastes.”
A special case of the derivation in J.M. Cushing’s text for the organism E. Coli2 is the set of nonlineardifferential equations
x′ =−0.075x+ (0.075)(0.005)− 1 63g(x)y, y′ =−0.075y+g(x)y,
(2)
where g(x) = 0.68x(0.0016 +x)−1. Of special interest to the study of this equation are two linearized equations at equilibria, given by
u′1 =−0.075u1−0.008177008175u2, u′2 = 0.4401515152u2,
(3)
v1′ =−1.690372243v1 −0.001190476190v2, v2′ = 101.7684513v1.
(4)
Although we cannot solve the nonlinear system explicitly, nevertheless there are explicit formulae foru1,u2,v1,v2 that complete the picture of how solutions x(t), y(t) behave at t= ∞. The result of the analysis is that E. Colisurvives indefinitely in this vessel at concentrationy ≈0.3.
Irregular Heartbeats and Lidocaine
The human malady of ventricular arrhythmiaor irregular heartbeat is treated clinically using the drug lidocaine.
Figure 7. Xylocaine label, a brand name for the drug lidocaine.
To be effective, the drug has to be maintained at a bloodstream concen- tration of 1.5 milligrams per liter, but concentrations above 6 milligrams per liter are considered lethal in some patients. The actual dosage de- pends upon body weight. The adult dosage maximum for ventricular tachycardia is reported at 3 mg/kg.3 The drug is supplied in 0.5%, 1%
and 2% solutions, which are stored at room temperature.
A differential equation model for the dynamics of the drug therapy uses
2In a biology Master’s thesis, two strains of Escherichia coli were grown in a glucose- limited chemostat coupled to a modified Robbins device containing plugs of silicone rubber urinary catheter material. Reference: Jennifer L. Adams and Robert J. C.
McLean, Applied and Environmental Microbiology, September 1999, p. 4285-4287, Vol. 65, No. 9.
3Source: Family Practice Notebook, http://www.fpnotebook.com/. The au- thor is Scott Moses, MD, who practises in Lino Lakes, Minnesota.
x(t) = amount oflidocaine in the bloodstream, y(t) = amount of lidocainein body tissue.
A typical set of equations, valid for a special body weight only, appears below; for more detail see J.M. Cushing’s text [?].
x′(t) =−0.09x(t) + 0.038y(t), y′(t) = 0.066x(t)−0.038y(t).
(5)
The physically significant initial data is zero drug in the bloodstream x(0) = 0 and injection dosage y(0) =y0. The answers:
x(t) =−0.3367y0e−0.1204t+ 0.3367y0e−0.0076t, y(t) = 0.2696y0e−0.1204t+ 0.7304y0e−0.0076t.
The answers can be used to estimate the maximum possible safe dosage y0 and the duration of time that the druglidocaine is effective.
Nutrient Flow in an Aquarium
Consider a vessel of water containing a radioactive isotope, to be used as a tracer for the food chain, which consists of aquatic plankton varieties Aand B.
Plankton are aquatic organisms that drift with the currents, typically in an environment like Chesapeake Bay. Plankton can be divided into two groups, phytoplankton and zooplankton. The phytoplankton are plant-likedrifters: diatoms and other alga. Zooplankton areanimal-like drifters: copepods, larvae, and small crustaceans.
Figure 8. Left: Bacillaria paxillifera, phytoplankton.
Right: Anomura Galathea zoea, zooplankton.
Let
x(t) = isotope concentration in the water, y(t) = isotope concentration in A,
z(t) = isotope concentration in B.
Typical differential equations are
x′(t) =−3x(t) + 6y(t) + 5z(t), y′(t) = 2x(t)−12y(t),
z′(t) =x(t) + 6y(t)−5z(t).
The answers are x(t) = 6c1+ (1 +√
6)c2e(−10+√6)t+ (1−√
6)c3e(−10−√6)t, y(t) =c1+c2e(−10+√6)t+c3e(−10−√6)t,
z(t) = 12
5 c1−2 +√
1.5c2e(−10+√6)t+−2 +√
1.5c3e(−10−√6)t. The constants c1, c2, c3 are related to the initial radioactive isotope concentrations x(0) = x0, y(0) = 0, z(0) = 0, by the 3×3 system of linear algebraic equations
6c1 + (1 +√
6)c2 + (1−√
6)c3 = x0,
c1 + c2 + c3 = 0,
12
5 c1 − 2 +√
1.5c2 + −2 +√
1.5c3 = 0.
Biomass Transfer
Consider a European forest having one or two varieties of trees. We select some of the oldest trees, those expected to die off in the next few years, then follow the cycle of living trees into dead trees. The dead trees eventually decay and fall from seasonal and biological events. Finally, the fallen trees become humus. Let variables x,y,z,t be defined by
x(t) = biomass decayed into humus, y(t) = biomass of dead trees,
z(t) = biomass of living trees,
t= time in decades (decade= 10 years).
A typical biological model is
x′(t) =−x(t) + 3y(t), y′(t) =−3y(t) + 5z(t), z′(t) =−5z(t).
Suppose there are no dead trees and no humus att= 0, with initiallyz0 units of living tree biomass. These assumptions imply initial conditions x(0) =y(0) = 0, z(0) =z0. The solution is
x(t) = 15
8 z0e−5t−2e−3t+e−t, y(t) = 5
2z0−e−5t+e−3t, z(t) =z0e−5t.
The live tree biomass z(t) =z0e−5t decreases according to a Malthusian decay law from its initial sizez0. It decays to 60% of its original biomass
in one year. Interesting calculations that can be made from the other formulae include the future dates when the dead tree biomass and the humus biomass are maximum. The predicted dates are approximately 2.5 and 8 years hence, respectively.
The predictions made by this model are trends extrapolated from rate observations in the forest. Like weather prediction, it is a calculated guess that disappoints on a given day and from the outset has no pre- dictable answer.
Total biomass is considered an important parameter to assess atmo- spheric carbon that is harvested by trees. Biomass estimates for forests since 1980 have been made by satellite remote sensing data with instances of 90% accuracy (Science87(5), September 2004).
Pesticides in Soil and Trees
A Washington cherry orchard was sprayed with pesticides.
Figure 9. Cherries in June.
Assume that a negligible amount of pesticide was sprayed on the soil.
Pesticide applied to the trees has a certain outflow rate to the soil, and conversely, pesticide in the soil has a certain uptake rate into the trees.
Repeated applications of the pesticide are required to control the insects, which implies the pesticide levels in the trees varies with time. Quantize the pesticide spraying as follows.
x(t) = amount of pesticide in the trees, y(t) = amount of pesticide in the soil,
r(t) = amount of pesticide applied to the trees, t= time in years.
A typical model is obtained from input-output analysis, similar to the brine tank models:
x′(t) = 2x(t)−y(t) +r(t), y′(t) = 2x(t)−3y(t).
In a pristine orchard, the initial data isx(0) = 0, y(0) = 0, because the trees and the soil initially harbor no pesticide. The solution of the model
obviously depends on r(t). The nonhomogeneous dependence is treated by the method of variation of parameters infra. Approximate formulae are
x(t)≈ Z t
0
1.10e1.6(t−u)−0.12e−2.6(t−u)r(u)du, y(t)≈
Z t 0
0.49e1.6(t−u)−0.49e−2.6(t−u)r(u)du.
The exponential rates 1.6 and −2.6 represent respectively the accumu- lation of the pesticide into the soil and the decay of the pesticide from the trees. The application rate r(t) is typically a step function equal to a positive constant over a small interval of time and zero elsewhere, or a sum of such functions, representing periodic applications of pesticide.
Forecasting Prices
A cosmetics manufacturer has a marketing policy based upon the price x(t) of its salon shampoo.
Figure 10. Salon shampoo sample.
The marketing strategy for the shampoo is to set the pricex(t) dynamically to reflect demand for the product. A low inventory is desirable, to reduce the overall cost of the product.
TheproductionP(t) and thesalesS(t) are given in terms of theprice x(t) and the change in price x′(t) by the equations
P(t) = 4− 3
4x(t)−8x′(t) (Production), S(t) = 15−4x(t)−2x′(t) (Sales).
The differential equations for the pricex(t) and inventory levelI(t) are x′(t) =k(I(t)−I0),
I′(t) =P(t)−S(t).
The inventory level I0 = 50 represents the desired level. The equations can be written in terms of x(t),I(t) as follows.
x′(t) = kI(t) − kI0,
I′(t) = 13
4 x(t) − 6kI(t) + 6kI0−11.
If k= 1, x(0) = 10 andI(0) = 7, then the solution is given by x(t) = 44
13 +86
13e−13t/2, I(t) = 50−43e−13t/2.
Theforecast of price x(t) ≈3.39 dollars at inventory levelI(t)≈50 is based upon the two limits
tlim→∞x(t) = 44
13, lim
t→∞I(t) = 50.
Coupled Spring-Mass Systems
Three masses are attached to each other by four springs as in Figure 11.
m1 m3
k2 k3 k4 k1
m2
Figure 11. Three masses
connected by springs. The masses slide along a frictionless horizontal surface.
The analysis uses the following constants, variables and assumptions.
Mass Constants
The massesm1,m2,m3are assumed to be point masses concentrated at their center of gravity.
Spring Constants
The mass of each spring is negligible. The springs op- erate according to Hooke’s law: Force = k(elongation).
Constantsk1,k2,k3,k4 denote the Hooke’s constants.
The springs restore after compression and extension.
Position Variables
The symbols x1(t), x2(t),x3(t) denote the mass posi- tions along the horizontal surface, measured from their equilibrium positions, plus right and minus left.
Fixed Ends The first and last spring are attached to fixed walls.
The competition method is used to derive the equations of motion.
In this case, the law is
Newton’s Second Law Force = Sum of the Hooke’s Forces.
The model equations are
m1x′′1(t) = −k1x1(t) +k2[x2(t)−x1(t)],
m2x′′2(t) = −k2[x2(t)−x1(t)] +k3[x3(t)−x2(t)], m3x′′3(t) = −k3[x3(t)−x2(t)]−k4x3(t).
(6)
The equations are justified in the case of all positive variables by observ- ing that the first three springs are elongated by x1, x2 −x1, x3 −x2, respectively. The last spring is compressed byx3, which accounts for the minus sign.
Another way to justify the equations is through mirror-image symmetry:
interchangek1 ←→k4,k2 ←→k3,x1←→x3, then equation 2 should be unchanged and equation 3 should become equation 1.
Matrix Formulation. System (6) can be written as a second order vector-matrix system
m1 0 0 0 m2 0 0 0 m3
x′′1 x′′2 x′′3
=
−k1−k2 k2 0 k2 −k2−k3 k3
0 k3 −k3−k4
x1 x2
x3
. More succinctly, the system is written as
Mx′′(t) =Kx(t)
where thedisplacement x,mass matrixM andstiffness matrixK are defined by the formulae
x=
x1 x2
x3
, M=
m1 0 0 0 m2 0 0 0 m3
, K=
−k1−k2 k2 0 k2 −k2−k3 k3
0 k3 −k3−k4
. Numerical example. Let m1 = 1, m2 = 1, m3 = 1, k1 = 2, k2 = 1, k3 = 1,k4 = 2. Then the system is given by
x′′1 x′′2 x′′3
=
−3 1 0 1 −2 1 0 1 −3
x1 x2
x3
. The vector solution is given by the formula
x1 x2 x3
= (a1cost+b1sint)
1 2 1
+a2cos√
3t+b2sin√ 3t
1 0
−1
+ (a3cos 2t+b3sin 2t)
1
−1 1
, where a1,a2,a3,b1,b2,b3 are arbitrary constants.
Boxcars
A special case of the coupled spring-mass system is three boxcars on a level track connected by springs, as in Figure 12.
k k
m m
m
Figure 12. Three identical boxcars connected by identical springs.
Except for the springs on fixed ends, this problem is the same as the one of the preceding illustration, therefore taking k1 =k4 = 0, k2=k3=k, m1 =m2 =m3=m gives the system
m 0 0 0 m 0 0 0 m
x′′1 x′′2 x′′3
=
−k k 0 k −2k k 0 k −k
x1
x2 x3
. Takek/m= 1 to obtain the illustration
x′′=
−1 1 0 1 −2 1 0 1 −1
x,
which has vector solution x = (a1+b1t)
1 1 1
+ (a2cost+b2sint)
1 0
−1
+a3cos√
3t+b3sin√ 3t
1
−2 1
, wherea1,a2,a3,b1,b2,b3 are arbitrary constants.
The solution expression can be used to discover what happens to the boxcars when the springs act normally upon compression but disengage upon expansion. An interesting physical situation is when one car moves along the track, contacts two stationary cars, then transfers its momen- tum to the other cars, followed by disengagement.
Electrical Network I
Consider theLR-network of Figure 13.
R1
i3 R3
R2
L3 L2
L1 i1
i2
Figure 13. An electrical network.
There are three resistorsR1,R2,R3
and three inductors L1,L2, L3. The currentsi1,i2,i3 are defined between nodes (black dots).
The derivation of the differential equations for the loop currents i1, i2, i3 uses Kirchhoff’s laws and the voltage drop formulae for resistors and inductors. The black dots in the diagram are thenodesthat determine the beginning and end of each of the currents i1, i2, i3. Currents are
defined only on the outer boundary of the network. Kirchhoff’s node law determines the currents across L2, L3 (arrowhead right) as i2−i1 and i3−i1, respectively. Similarly,i2−i3 is the current acrossR1 (arrowhead down). Using Ohm’s law VR = RI and Faraday’s law VL = LI′ plus Kirchhoff’s loop law algebraic sum of the voltage drops is zero around a closed loop (see themaple code below), we arrive at the model
i′1 = −
R2 L1
i2 −
R3 L1
i3, i′2 = −
R2 L2 +R2
L1
i2 +
R1 L2 − R3
L1
i3, i′3 =
R1 L3 −R2
L1
i2 − R1
L3 +R3 L1 +R3
L3
i3
A computer algebra system is helpful to obtain the differential equations from the closed loop formulas. Part of the theory is that the number of equations equals the number ofholesin the network, called theconnec- tivity. Here’s somemaple code for determining the equations in scalar and also in vector-matrix form.
loop1:=L1*D(i1)+R3*i3+R2*i2=0;
loop2:=L2*D(i2)-L2*D(i1)+R1*(i2-i3)+R2*i2=0;
loop3:=L3*D(i3)-L3*D(i1)+R3*i3+R1*(i3-i2)=0;
f1:=solve(loop1,D(i1));
f2:=solve(subs(D(i1)=f1,loop2),D(i2));
f3:=solve(subs(D(i1)=f1,loop3),D(i3));
with(linalg):
jacobian([f1,f2,f3],[i1,i2,i3]);
Electrical Network II
Consider theLR-network of Figure 14. This network produces only two differential equations, even though there are three holes (connectivity 3). The derivation of the differential equations parallels the previous network, so nothing will be repeated here.
A computer algebra system is used to obtain the differential equations from the closed loop formulae. Below is maple code to generate the equations i′1 =f1,i′2=f2,i3=f3.
loop1:=L1*D(i1)+R2*(i1-i2)+R1*(i1-i3)=0;
loop2:=L2*D(i2)+R3*(i2-i3)+R2*(i2-i1)=0;
loop3:=R3*(i3-i2)+R1*(i3-i1)=E;
f3:=solve(loop3,i3);
f1:=solve(subs(i3=f3,loop1),D(i1));
f2:=solve(subs(i3=f3,loop2),D(i2));
E
R1 R2
i1 L1
R3
i3 i2
L2
Figure 14. An electrical network.
There are three resistorsR1,R2,R3, two inductorsL1,L2 and a batteryE.
The currentsi1,i2,i3 are defined between nodes (black dots).
The model, in the special case L1=L2= 1 and R1=R2=R3 =R:
i′1 = − 3R
2 i1 + 3R
2 i2 + E 2, i′2 = 3R
2 i1 − 3R
2 i2 + E 2,
i3 = 1
2i1 + 1
2i2 + E 2R.
It is easily justified that the solution of the differential equations for initial conditionsi1(0) =i2(0) = 0 is given by
i1(t) = E
2t, i2(t) = E 2t.
Logging Timber by Helicopter
Certain sections of National Forest in the USA do not have logging ac- cess roads. In order to log the timber in these areas, helicopters are employed to move the felled trees to a nearby loading area, where they are transported by truck to the mill. The felled trees are slung beneath the helicopter on cables.
Figure 15. Helicopter logging.
Left: An Erickson helicopter lifts felled trees.
Right: Two trees are attached to the cable to lower transportation costs.
The payload for two trees approximates a double pendulum, which oscil- lates during flight. The angles of oscillationθ1,θ2 of the two connecting cables, measured from the gravity vector direction, satisfy the following differential equations, in which g is the gravitation constant, m1, m2
denote the masses of the two trees and L1,L2 are the cable lengths.
(m1+m2)L21θ′′1 + m2L1L2θ2′′ + (m1 +m2)L1gθ1= 0, m2L1L2θ′′1 + m2L22θ2′′ + m2L2gθ2= 0.
This model is derived assuming small displacements θ1, θ2, that is, sinθ≈θ for both angles, using the following diagram.
θ2
L1 L2 m2
m1 θ1
Figure 16. Logging Timber by Helicopter.
The cables have lengthsL1, L2. The anglesθ1,θ2are measured from vertical.
The lengths L1,L2 are adjusted on each trip for the length of the trees, so that the trees do not collide in flight with each other nor with the helicopter. Sometimes, three or more smaller trees are bundled together in a package, which is treated here as identical to a single, very thick tree hanging on the cable.
Vector-matrix model. The angles θ1, θ2 satisfy the second-order vector-matrix equation
(m1+m2)L1 m2L2
L1 L2
! θ1 θ2
!′′
=− m1g+m2g 0
0 g
! θ1 θ2
! . This system is equivalent to the second-order system
θ1 θ2
!′′
=
−m1g+m2g L1m1
m2g L1m1 m1g+m2g
L2m1 −(m1+m2)g L2m1
θ1 θ2
! .
Earthquake Effects on Buildings
A horizontal earthquake oscillation F(t) =F0cosωtaffects each floor of a 5-floor building; see Figure 17. The effect of the earthquake depends upon the natural frequencies of oscillation of the floors.
In the case of a single-floor building, the center-of-mass position x(t) of the building satisfies mx′′+kx = E and the natural frequency of oscillation ispk/m. The earthquake forceEis given by Newton’s second law: E(t) =−mF′′(t). Ifω≈pk/m, then the amplitude ofx(t) is large compared to the amplitude of the force E. The amplitude increase in x(t) means that a small-amplitude earthquake wave can resonant with the building and possibly demolish the structure.
3
F 4 5
1 2
Figure 17. A 5-Floor Building.
A horizontal earthquake waveF affects every floor. The actual wave has wavelength many times larger than the illustration.
The following assumptions and symbols are used to quantize the oscilla- tion of the 5-floor building.
• Each floor is considered a point mass located at its center-of-mass.
The floors have masses m1, . . . , m5.
• Each floor is restored to its equilibrium position by a linear restor- ing force or Hooke’s force −k(elongation). The Hooke’s constants are k1, . . . , k5.
• The locations of masses representing the 5 floors are x1, . . . , x5. The equilibrium position is x1 =· · ·=x5 = 0.
• Damping effects of the floors are ignored. This is a frictionless system.
The differential equations for the model are obtained bycompetition:
the Newton’s second law force is set equal to the sum of the Hooke’s forces and the external force due to the earthquake wave. This results in the following system, wherek6 = 0, Ej =mjF′′ for j = 1,2,3,4,5 and F =F0cosωt.
m1x′′1 = −(k1+k2)x1+k2x2+E1, m2x′′2 = k2x1−(k2+k3)x2+k3x3+E2, m3x′′3 = k3x2−(k3+k4)x3+k4x4+E3, m4x′′4 = k4x3−(k4+k5)x4+k5x5+E4, m5x′′5 = k5x4−(k5+k6)x5+E5.
In particular, the equations for a floor depend only upon the neighboring floors. The bottom floor and the top floor are exceptions: they have just one neighboring floor.
Vector-matrix second order system. Define
M =
m1 0 0 0 0
0 m2 0 0 0
0 0 m3 0 0
0 0 0 m4 0
0 0 0 0 m5
, x=
x1 x2
x3 x4 x5
, H=
E1 E2
E3 E4 E5
,
K=
−k1−k2 k2 0 0 0
k2 −k2−k3 k3 0 0
0 k3 −k3−k4 k4 0
0 0 k4 −k4−k5 k5
0 0 0 k5 −k5−k6
.
In the last row, k6 = 0, to reflect the absence of a floor above the fifth.
The second order system is
Mx′′(t) =Kx(t) +H(t).
The matrixM is called themass matrixand the matrixK is called the Hooke’s matrix. The external force H(t) can be written as a scalar functionE(t) =−F′′(t) times a constant vector:
H(t) =−ω2F0cosωt
m1
m2 m3 m4
m5
.
Identical floors. Let us assume that all floors have the same mass m and the same Hooke’s constant k. ThenM =mI and the equation becomes
x′′=m−1
−2k k 0 0 0
k −2k k 0 0
0 k −2k k 0
0 0 k −2k k
0 0 0 k −k
x−F0ω2cos(ωt)
1 1 1 1 1
.
The Hooke’s matrix K is symmetric (KT = K) with negative entries only on the diagonal. The last diagonal entry is −k (a common error is to write −2k).
Particular solution. The method of undetermined coefficients predicts a trial solution xp(t) = c cosωt, because each differential equation has nonhomogeneous term −F0ω2cosωt. The constant vector c is found by trial solution substitution. Cancel the common factor cosωt in the substituted equation to obtain the equation m−1K+ω2Ic=F0ω2b, where b is column vector of ones in the preceding display. Let B(ω) = m−1K+ω2I. Then the formula B−1= adj(B)
det(B) gives c=F0ω2adj(B(ω))
det(B(ω))b.
The constant vector ccan have a large magnitude when det(B(ω))≈0.
This occurs when −ω2 is nearly an eigenvalue of m−1K.
Homogeneous solution. The theory of this chapter gives the homo- geneous solutionxh(t) as the sum
xh(t) =
5
X
j=1
(ajcosωjt+bjsinωjt)vj wherer=ωj and v=vj 6=0 satisfy
1
mK+r2I
v=0.
Special case k/m= 10. Then
1 mK =
−20 10 0 0 0
10 −20 10 0 0
0 10 −20 10 0
0 0 10 −20 10
0 0 0 10 −10
and the valuesω1, . . . ,ω5 are found by solving the determinant equation det((1/m)K+ω2I) = 0, to obtain the values in Table 1.
Table 1. The natural frequencies for the special casek/m= 10.
Frequency Value
ω1 0.900078068
ω2 2.627315231
ω3 4.141702938
ω4 5.320554507
ω5 6.068366391
General solution. Superposition implies x(t) = xh(t) +xp(t). Both terms of the general solution represent bounded oscillations.
Resonance effects. The special solution xp(t) can be used to obtain some insight into practical resonance effects between the incoming earth- quake wave and the building floors. When ω is close to one of the fre- quencies ω1, . . . , ω5, then the amplitude of a component of xp can be very large, causing the floor to take an excursion that is too large to maintain the structural integrity of the floor.
Thephysical interpretationis that an earthquake wave of the proper frequency, having time duration sufficiently long, can demolish a floor and hence demolish the entire building. The amplitude of the earthquake wave does not have to be large: a fraction of a centimeter might be enough to start the oscillation of the floors.
