Issuu

SOLUTION MANUAL

Part 1 - Modeling, Computers, and Error Analysis Chapter

1) Mathematical Modeling and Engineering Problem Solving

Chapter

2) Programming and Software

Chapter

3) Approximations and Round-Off Errors

Chapter

4) Truncation Errors and the Taylor Series

Part 2 - Roots of Equations Chapter

5) Bracketing Methods

Chapter

6) Open Methods

Chapter

7) Roots of Polynomials

Chapter

8) Case Studies: Roots of Equations

Part 3 - Linear Algebraic Equations Chapter

9) Gauss Elimination

Chapter

10) LU Decomposition and Matrix Inversion

Chapter

11) Special Matrices and Gauss-Seidel

Chapter

12) Case Studies: Linear Algebraic Equations

Part 4 - Optimization Chapter

13) One-Dimensional Unconstrained Optimization

Chapter

14) Multidimensional Unconstrained Optimization

Chapter

15) Constrained Optimization

Chapter

16) Case Studies: Optimization

Part 5 - Curve Fitting Chapter

17) Least-Squares Regression

Chapter

18) Interpolation

Chapter

19) Fourier Approximation

Chapter

20) Case Studies: Curve Fitting

Part 6 - Numerical Differentiation and Integration

Chapter

21) Newton-Cotes Integration Formulas

Chapter

22) Integration of Equations

Chapter

23) Numerical Differentiation

Chapter

24) Case Studies: Numerical Integration and Differentiation

Part 7 - Ordinary Differential Equations Chapter

25) Runge-KuttaMethods

Chapter

26) Stiffness and Multistep Methods

Chapter

27) Boundary-Value and Eigenvalue Problems

Chapter

28) Case Studies: Ordinary Differential Equations

Part 8 - Partial Differential Equations Chapter

29) Finite Difference: Elliptic Equations

Chapter

30) Finite Difference: Parabolic Equations

Chapter

31) Finite-Element Method

Chapter

32) Case Studies: Partial Differential Equations.

CHAPTER 1 Use calculus to solve Eq. (1.9) for the case where the initial velocity υ(0) is nonzero.

1.1

We will illustrate two different methods for solving this problem: (1) separation of variables, and (2) Laplace transform. dv

 g

Separation of variables: Separation of variables gives



1 c

g

dv  v

 dt

The integrals can be evaluated as



c   ln  g  v  m  

tC

c/m

where C = a constant of integration, which can be evaluated by applying the initial condition to yield c   ln  g  v (0)  m  

C 

c/m

which can be substituted back into the solution



c   ln  g  v  m  

t

c   ln  g  v (0)  m  

c/m

This result can be rearranged algebraically to solve for v, v  v (0) e

 (c / m )t



mg c

1  e

 (c / m )t



where the first part is the general solution and the second part is the particular solution for the constant forcing function due to gravity. For the case where, v(0) = 0, the solution reduces to Eq. (1.10) v

mg c

1  e

 (c / m )t



Laplace transform solution: An alternative solution is provided by applying Laplace transform to the differential equation to give sV ( s )  v (0) 

g s



V (s)

Solve algebraically for the transformed velocity V (s) 

v (0) sc/m



(1)

s(s  c / m )

The second term on the right of the equal sign can be expanded with partial fractions g s(s  c / m)





sc/m

A ( s  c / m )  Bs



s(s  c / m )

(2)

By equating like terms in the numerator, the following must hold g  A

0  As  Bs

The first equation can be solved for A = mg/c. According to the second equation, B = –A, so B = –mg/c. Substituting these back into (2) gives g s(s  c / m )



mg / c



mg / c sc/m

This can be substituted into Eq. 1 to give V (s) 

v (0) sc/m



mg / c



mg / c sc/m

Taking inverse Laplace transforms yields v ( t )  v (0) e

 (c / m )t





1  e

 (c / m )t

or collecting terms v ( t )  v (0) e

 (c / m )t



 (c / m )t



1.2 Repeat Example 1.2. Compute the velocity to t = 10 s, with a step size of (a) 1 and (b) 0.5 s. Can you make any statement regarding the errors of the calculation based on the results? At t = 10 s, the analytical solution is 44.91893 (Example 1.1). The relative error can be calculated with relative true error 

analytical  num erical

 100%

analytical

The numerical results are: step

v(10)

2 1 0.5

48.0179 46.4112 45.6509

magnitude of relative error 6.899% 3.322% 1.630%

The error versus step size can then be plotted as Copyright 2021 © McGraw-Hill Education. All rights reserved. No reproduction or distribution without the prior written consent of McGraw-Hill Education.

Relative True Error vs. Step size 8.0 Relative error (%)

7.0 6.0 5.0 4.0 3.0 2.0 1.0 0.0 0

0.5

1.5

2.5

Step size (seconds)

Thus, halving the step size approximately halves the error. 1.3 Rather than the linear relationship of Eq. (1.7), you might choose to model the upward force on the parachutist as a second-order relationship, FU = - c ¢v ∣ v ∣

where c′ = a bulk second-order drag coefficient (kg/m). Note that the second-order term could be represented as v2 if the parachutist always fell in the downward direction. For the present case, we use the more general representation, v ∣ v ∣ , so that the proper sign is obtained for both the downward and the upward directions. (a) Using calculus, obtain the closed-form solution for the case where the jumper is initially at rest (υ = 0 at t = 0). (b) Repeat the numerical calculation in Example 1.2 with the same initial condition and parameter values, but with second-order drag. Use a value of 0.225 kg/m for cd′. (a) You are given the following differential equation with the initial condition, v(t = 0) = 0, dv

 g

c

Multiply both sides by m/c′ gives m dv c  dt



m c

g v

Define a 

mg / c

m dv

c  dt

 a v

Integrate by separation of variables, c

 a  v   m dt 2

A table of integrals can be consulted to find that dx

 a  x  a tanh 2

1 x

Therefore, the integration yields 1

tanh

1 v



c

tC

If v = 0 at t = 0, then because tanh–1(0) = 0, the constant of integration C = 0 and we obtain the equation 1

tanh

1 v



c

This result can then be rearranged to solve for v  tanh   c 

v

gc   t m 

(b) Using Euler’s method, the first two steps are computed 0.225  2 v (2)  0   9.81  (0)  2  19.62 68.1   0.225  2 v (4)  19.62   9.81  (1 9.62)  2  3 6.696 3 1 4 5 4 6 8. 1  

The computation can be continued and the results summarized along with the analytical result as: t

v-numerical

0 2 4 6 8 10 12 

0 19.62 36.69631454 47.41794779 52.18027088 53.80832813 54.29615076 54.48999908

dv/dt 9.81 8.538157 5.360817 2.381162 0.814029 0.243911 0.069674 0

v-analytical 0 18.8138836 33.61984724 43.22542283 48.7004867 51.59332241 53.06072073 54.48999908

A plot of the numerical and analytical results can be developed

Second Order Drag 60

velocity (m/s)

50 40 30

Euler Approx dt = 2

Analytical

10 0 0

time (s)

1.4 For the free-falling parachutist with linear drag, assume a first jumper is 70 kg and has a drag coefficient of 12 kg/s. If a second jumper has a drag coefficient of 15 kg/s and a mass of 75 kg, how long will it take him to reach the same velocity the first jumper reached in 10 s? v (t ) 

(1  e

 (c / m ) t

)

Solve the equation for time as a function of velocity t= -

æ ö vc ÷ ÷ ln çç1 ÷ çè c mg ÷ ø

jumper #1: v ( t ) 

9.81(70)

(1  e

 (12 / 70 ) 9

)  4 6 .91922

jumper #2: 46.91922  t= -

9.81(75)

(1  e

 (15 / 75 ) t

)

æ 15(46.91922) ö ÷ ÷ ln çç1 ÷ = 15.68175 15 çè 75(9.81) ÷ ø 75

The second jumper will reach the t=10 s velocity of the first jumper after about 15.68 seconds. 1.5 Compute the velocity of a parachutist using Euler’s method for the case where m = 80 kg and c = 10 kg/s. Perform the calculation from t = 0 to 20 s with a step size of 1 s. Use an initial condition that the parachutist has an upward velocity of 20 m/s at t = 0. At t = 10 s, assume that the chute is instantaneously deployed so that the drag coefficient jumps to 50 kg/s. Before the chute opens (t < 10), Euler’s method can be implemented as 10   v ( t   t )  v ( t )   9.81  v (t )   t 80  

6 After the chute opens (t  10), the drag coefficient is changed and the implementation becomes 50   v ( t   t )  v ( t )   9.81  v (t )   t 80  

Here is a summary of the results along with a plot:

t 0 1 2 3 4 5 6 7 8 9

Chute closed dv/dt v -20.00000 12.31000 -7.69000 10.77125 3.08125 9.42484 12.50609 8.24674 20.75283 7.21590 27.96873 6.31391 34.28264 5.52467 39.80731 4.83409 44.64139 4.22983 48.87122 3.70110

t 10 11 12 13 14 15 16 17 18 19 20

Chute opened dv/dt v 52.57232 -23.04770 29.52462 -8.64289 20.88173 -3.24108 17.64065 -1.21541 16.42524 -0.45578 15.96947 -0.17092 15.79855 -0.06409 15.73446 -0.02404 15.71042 -0.00901 15.70141 -0.00338 15.69803 -0.00127

Parachutist velocity from rising balloon

velocity (m/s)

60 50

Euler approximation dt=1 s

Analytical Value

30 20 10 0 -10 -20

Chute deploys at t=10

-30

1.6

time (s)

The following information is available for a bank account: Date

Deposits

Withdrawals

220.13

327.26

216.80

378.61

450.25

106.80

127.31

350.61

5/1

Interest

Balance 1512.33

6/1 7/1 8/1 9/1

Note that the money earns interest, which is computed as Interest = iB i

where i = the interest rate expressed as a fraction per month and Bi = the initial balance at the beginning of the month. (a) Use the conservation of cash to compute the balance on 6/1, 7/1, 8/1, and 9/1 if the interest rate is 1% per month (i = 0.01/month). Show each step in the computation. (b) Write a differential equation for the cash balance in the form dB

= f ( D ( t ), W ( t ), i )

where t = time (months), D(t) = deposits as a function of time ($/month), W(t) = withdrawals as a function of time ($/month). For this case, assume that interest is compounded continuously; that is, interest = iB. (c) Use Euler’s method with a time step of 0.5 month to simulate the balance. Assume that the deposits and withdrawals are applied uniformly over the month. (d) Develop a plot of balance versus time for (a) and (c). (a) This is a transient computation. For the period ending June 1: Balance = Previous Balance + Deposits – Withdrawals + Interest Balance = 1522.33 + 220.13 – 327.26 + 0.01(1522.33) = 1430.42 Note that the interest added to the account is for maintaining the account the previous month; for example, the interest collected during from May 1 to June 1 is calculated on the May 1 st balance. The balances for the remainder of the periods can be computed in a similar fashion as tabulated below: Date 1-May

Deposit

Withdrawal

Interest

$ 220.13

$ 327.26

15.12

$ 216.80

$ 378.61

14.20

$ 450.25

$ 106.80

12.73

$ 127.31

$ 350.61

16.29

1-Jun

$ 1,420.32

1-Jul

$ 1,272.72

1-Aug

$ 1,628.89

1-Sep

(b)

Balance $ 1,512.33

$ 1,421.88

 D ( t )  W ( t )  iB

 220.13  327.26  0.01(1512.33)   92.01

dt B (0.5)  1512.33  92.01(0.5)  1466.33

for t = 0.5 to 1: dB

 220.13  327.260  0.01(1466.33)   92.47

dt B (0.5)  1466.33  92.47(0.5)  1420.09

The balances for the remainder of the periods can be computed in a similar fashion as tabulated below. Parenthesis indicate negative numbers. Date 1-May 16-May 1-Jun 16-Jun 1-Jul 16-Jul 1-Aug 16-Aug 1-Sep

Deposit

Withdrawal

$ 220.13 $ $ 220.13 $ $ 216.80 $ $ 216.80 $ $ 450.25 $ $ 450.25 $ $ 127.31 $ $ 127.31 $

Interest

327.26 $ 327.26 $ 378.61 $ 378.61 $ 106.80 $ 106.80 $ 350.61 $ 350.61 $

dB/dt

Balance

15.12 $ (92.01) $ 1,512.33 14.66 $ (92.47) $ 1,466.33 14.20 $ (147.61) $ 1,420.09 13.46 $ (148.35) $ 1,346.29 12.72 $ 356.17 $ 1,272.12 14.50 $ 357.95 $ 1,450.20 16.29 $ (207.01) $ 1,629.18 15.26 $ (208.04) $ 1,525.67 $ 1,421.65

(d) As in the plot below, the results of the two approaches are very close.

Compounded Interest $1,800.00 Bi-monthly Interest

$1,700.00

Monthly Interest

$1,600.00 $1,500.00 $1,400.00 $1,300.00 $1,200.00 5-Apr

25-Apr

15-May

4-Jun

24-Jun

14-Jul

3-Aug

23-Aug

12-Sep

1.7 The amount of a uniformly distributed radioactive contaminant contained in a closed reactor is measured by its concentration c (becquerel/liter, or Bq/L). The contaminant decreases at a decay rate proportional to its concentration—that is, decay rate = - kc

where k is a constant with units of day−1. Therefore, according to Eq. (1.13), a mass balance for the reactor can be written as dc

- kc

dt æchange ö ÷ = çç ÷ ÷ çèin m ass ø

ædecrease ö ÷ çç ÷ ÷ èç by decay ø

(a) Use Euler’s method to solve this equation from t = 0 to 1 d with k = 0.2 d−1. Employ a step size of Δt = 0.1. The concentration at t = 0 is 10 Bq/L. (b) Plot the solution on a semilog graph (i.e., ln c versus t) and determine the slope. Interpret your results.

(a) The first two steps are c (0.1)  100  0.2(10)0.1  9.8 Bq/L c (0.2)  9.80  0.2(9.80)0.1  9.604 Bq/L

The process can be continued to yield t

dc/dt -2.0000 -1.9600 -1.9208 -1.8824 -1.8447 -1.8078 -1.7717 -1.7363 -1.7015 -1.6675 -1.6341

c 10.0000 9.8000 9.6040 9.4119 9.2237 9.0392 8.8584 8.6813 8.5076 8.3375 8.1707

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1

(b) The results when plotted on a semi-log plot yields a straight line

Euler concentration 2.35 Euler concentration

ln ([c]) Bq/L

2.3

Linear (Euler concentration)

2.25 2.2 y = -0.202027x + 2.302585 R² = 1.000000

2.15 2.1 2.05 0

0.2

0.4

0.6

0.8

1.2

time (days)

The slope of this line can be estimated from the first and last points as ln(8.1707)  ln(10)

  0.2020

Thus, the slope is approximately equal to the negative of the decay rate, within 10%. If we had used a smaller step size, the result would be more exact. A linear regression fit is shown on the graph using all points.

1.8 A group of 35 students attend a class in a room that measures 11 m by 8 m by 3 m. Each student takes up about 0.075 m3 and gives out about 80 W of heat (1 W = 1 J/s). Calculate the air temperature rise during the first 15 minutes of the class if the room is completely sealed and insulated. Assume the heat capacity, Cυ, for air is 0.718 kJ/(kg K). Assume air is an ideal gas at 20°C and 101.325 kPa. Note that the heat absorbed by the air Q is related to the mass of the air m, the heat capacity, and the change in temperature by the following relationship: Q  m

 C dT  m C (T  T ) 2

The mass of air can be obtained from the ideal gas law: m

PV 

Mwt

where P is the gas pressure, V is the volume of the gas, Mwt is the molecular weight of the gas (for air, 28.97 kg/kmol), and R is the ideal gas constant [8.314 kPa m3/(kmol K)]. Q students  35 ind  80 m 

PV M w t

J ind s

s m in



 2, 520 kJ

1000 J 3



(101.325 kPa )(11m  8m  3m  35  0.075 m )(28.97 kg/km ol) 3

(8.314 kPa m / (km ol K )((20  273.15)K )

T 

 15 m in  60

Q students

2, 520 kJ



mCv

 314.796 kg

 11.14928 K

(314.796 kg )(0.718 kJ/(kg K ))

The final temperature is 20 + 11.14928 = 31.14928 oC, though only  T is requested. 1.9 A storage tank contains a liquid at depth y, where y = 0 when the tank is half full (Fig. P1.9). Liquid is withdrawn at a constant flow rate Q to meet demands. The contents are resupplied at a sinusoidal rate 3Q sin2(t).

Equation (1.13) can be written for this system as d ( Ay)

3 Q sin ( t )

æchange in ö ÷ çç = ÷ ÷ èç volum e ø

(inflow )

(outflow )

= 3

sin ( t ) -

Q A

Use Euler’s method to solve for the depth y from t = 0 to 10 d with a step size of 0.5 d. The parameter values are A = 1200 m2 and Q = 500 m3/d. Assume that the initial condition is y = 0. The first two steps yield 500   500 2 y (0.5)  0   3 sin (0)   0.5  0  (  0.41667 ) 0.5   0.20833 1200 1200   500   500 2 y (1)   0.20833   3 sin (0.5)  0.5   0.20833  (  0.12936) 0.5   0.23588 1200   1200

The process can be continued to give the following table and plot: t 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5

y 0.00000 -0.20833 -0.27301 -0.03880 0.37474 0.68317 0.69869 0.50281 0.37138 0.52101 0.90991

dy/dt -0.41667 -0.12936 0.46843 0.82708 0.61686 0.03104 -0.39177 -0.26286 0.29927 0.77779 0.73275

t 5.5 6 6.5 7 7.5 8 8.5 9 9.5 10

dy/dt 0.20557 -0.31908 -0.35882 0.12287 0.68314 0.80687 0.38031 -0.20436 -0.40961 -0.04672

y 1.27629 1.37907 1.21953 1.04012 1.10156 1.44313 1.84656 2.03672 1.93453 1.72973

graph of result of Euler. Analytical solution also shown as a cross-check. Smaller step sizes match better.

sinusoidal inflow of a storage tank 2.5 Euler Approximation

water height (m)

2.0

Analytical Solution

1.5 1.0 0.5 0.0 0 -0.5

time(days)

1.10 For the same storage tank described in Prob. 1.9, suppose that the outflow is not constant but rather depends on the depth. For this case, the differential equation for depth can be written as dy

= 3

sin ( t ) -

a (1 + y )

1.5

Use Euler’s method to solve for the depth y from t = 0 to 10 d with a step size of 0.5 d. The parameter values are A = 1200 m2, Q = 500 m3/d, and α = 150. Assume that the initial condition is y = 0. The first two steps yield 1.5  500 150(1  0)  2 y (0.5)  0   3 sin (0)   0.5  0  0.12(0.5)   0.0625 1200  1200  1.5  500 150(1  0.0625)  2 y (1)   0.0625   3 sin (0.5)   0.5   0.0625  0.17384(0.5)  0 .02442 1200  1200 

The process can be continued to give t 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5

y 0.00000 -0.06250 0.02442 0.40216 0.92027 1.27072 1.28072 1.07789 0.96759 1.15306 1.55284

dy/dt -0.12500 0.17384 0.75548 1.03620 0.70090 0.02000 -0.40565 -0.22060 0.37094 0.79955 0.63957

t 5.5 6 6.5 7 7.5 8 8.5 9 9.5 10

dy/dt 0.01364 -0.51317 -0.47313 0.07876 0.62757 0.65677 0.12530 -0.48005 -0.60721 -0.15013

y 1.87262 1.87944 1.62286 1.38629 1.42567 1.73946 2.06784 2.13049 1.89046 1.58686

time-dependent inflow + outflow 2.0

water height (m)

1.5 1.0 Euler Aprroximation

0.5 0.0 0 -0.5

time(days)