sm ch (21).pdf

1

CHAPTER 21 21.1 x xi2 xi1 xi xi+1 xi+2

f ( x) x)

0.261799388 0.523598776 0.785398163 1.047197551 1.308996939

0.965925826 0.866025404 0.707106781 0.5 0.258819045

true = sin( sin( /4) /4) = 0.70710678 The results are summarized as first-order Forward Backward Centered

0.79108963 11.877% 0.60702442 14.154% 0.69905703 1.138%

second-order 0.72601275 2.674% 0.71974088 1.787% 0.70699696 0.016%

21.2 x xi2 xi1 xi xi+1 xi+2

f ( x) x)

1.8 1.9 2 2.1 2.2

6.049647464 6.685894442 7.389056099 8.166169913 9.025013499

Both the first and second derivatives have the same value, truth  e 2  7.389056099 The results are summarized as

First derivative Second derivative

first-order 7.401377351 -0.166750% 7.395215699 -0.083361%

second-order 7.389031439 0.000334% 7.389047882 0.000111%

21.3 First, we will use forward expansions. The Taylor series expansion about a = xi and x = xi+2 (2 x steps forward) can be written as: f ( xi  2 )  f ( xi )  f '( xi )2x 

1 2

f "( xi )(2x) 2 

f ( xi  2 )  f ( xi )  2 f '( xi ) x  2 f "( xi ) x 2 

8 6

1 6

f ( xi )(2x)3 

f ( xi )x 3 

16 24

1 24

f

( 4)

( xi )(2 x) 4 

1 120

f (5) ( xi )(2 x) 5    

f ( 4) ( xi )x 4



32 120

f

( 5)

(1)

5

( xi ) x    

PROPRIETARY MATERIAL. © The McGraw-Hill Companies, Inc. All rights rights reserved. No part of this Manual may be displayed, reproduced or distributed in any form or by any means, without the prior written permission of the publisher, or used beyond the limited distribution to teachers and educators permitted by McGraw-Hill for their individual course preparation. If you are a student using this Manual, you are using it without permission. permission.

2 Taylor series expansion about a = xi and x = xi+1 ( x steps forward): f ( xi 1 )  f ( xi )  f '( '( xi )x 

1 2

f "( xi )x 2 

1

f ( xi ) x 3 

6

1

f

24

(4 )



( xi ) x 4

1 120

f

(5)

(2)

(3)

5

( xi ) x    

Multiply Eq. 2 by 2 and subtract the result from Eq. 1 to yield f ( xi  2 )  2 f ( xi 1 )   f ( xi )  f "( "( xi ) x 2 

6

f ( xi ) x 3 

6

14

f ( 4) ( x i ) x 4

24



30 120

f

(5)

5

( xi ) x    

Next, we will use backward expansions. The Taylor series expansion about a = xi and x = xi2 (2  x steps backward) can be written as: f ( xi  2 )  f ( xi )  f '( xi )( 2x) 

1 2

f "( xi )( 2 x) 2 

1 6

f ( xi )( 2 x) 3 



f ( xi  2 )  f ( xi )  2 f '( xi ) x  2 f "( xi )x 2 

8 6

16

f ( xi ) x 3 

24

1 120

1 24

f

( 4)

( xi )( 2 x) 4

f (5) ( xi )(  2x)5 

f ( 4) ( xi ) x 4

-

32 120

f

(5)

(4)

5

( xi ) x    

Taylor series expansion about a = xi and x = xi1 ( x steps backward): f ( xi 1 )  f ( xi )  f '( '( xi ) x 

1 2

f "( xi ) x 2 

1 6

f ( xi ) x 3 

1

f ( 4 ) ( xi ) x 4

24



1 120

f

( 5)

(5)

(6)

5

( xi ) x    

Multiply Eq. 5 by 2 and subtract the result from Eq. 4 to yield 2 f ( xi 1 )  f ( xi  2 )  f ( xi )  f "( "( xi ) x 2 

6 6

f ( xi ) x 3 

14 24

f

(4)

( x i ) x 4



30 120

f

(5)

5

( xi ) x    

Add Eqs (3) and (6) f ( xi  2 )  2 f ( xi 1 )  2 f ( xi 1 )  f ( xi 2 )  2 f ( xi ) x 3 

60 120

f

(5)

( xi ) x 5    

(7)

Equation 7 can be solved for

f ( xi ) 

f ( xi  2 )  2 f ( xi 1 )  2 f ( xi1 )  f ( xi 2 )

2 x

3

60

 120

f (5) ( xi )x5

2x 3

 

PROPRIETARY MATERIAL. © The McGraw-Hill Companies, Inc. All rights rights reserved. No part of this Manual may be displayed, reproduced or distributed in any form or by any means, without the prior written permission of the publisher, or used beyond the limited distribution to teachers and educators permitted by McGraw-Hill for their individual course preparation. If you are a student using this Manual, you are using it without permission. permission.

3

f ( xi ) 

f ( xi  2 )  2 f ( xi 1 )  2 f ( xi1 )  f ( xi 2 )

2 x3

1

 f (5) ( xi )x 2   4

21.4 The true value is sin( /4) = 0.70710678. D(  / 3)  D(  / 6)  D 

0.25882  0.965926

  0.58477

2(1.047198) 0.258819  0.965926 2(0.523599)

 0.67524

4

1 (0.67524)  (0.58477)  0.70539 3 3

21.5 The true value is 1/ x = 1/5 = 0.2.

1.94591  1.098612

D(2)  D(1)  D 

2(2)

 0.211824

1.791759  1.386294 2(1)

 0.202733

4

1 (0.202733)  (0.211824)  0.199702 3 3

21.6 The true value 3

2

f '(0)  8(0)  18(0)  12  12

Equation (21.21) can be used to compute the derivative as x0 = 0.5 x1 = 1 x2 = 2 f '(0)   1.125

f ( x0) = 1.125 f ( x1) = 24 f ( x2) = 48

2(0)  1  2 ( 0.5  1)(0.5  2)

 ( 24)

2(0)  ( 0.5)  2 (1  ( 0.5))(1  2)

 ( 48)

2(0)  ( 0.5)  1 (2  ( 0.5))(2  1)

 0.9  24  9.6  13.5

Centered difference: f '(0) 

24  12   18 1  (1)

21.7 At x = xi, Eq. (21.21) is f '( x )  f ( xi 1 )

2 xi  xi  xi 1 ( xi 1  xi )( xi1  xi1 )

 f ( xi )

2 xi  xi 1  xi1 ( xi  xi1 )( xi  xi1 )

 f ( xi 1 )

2 xi  xi1  xi ( xi 1  xi 1 )( x i1  xi )

For equispaced points that are h distance apart, this equation becomes

PROPRIETARY MATERIAL. © The McGraw-Hill Companies, Inc. All rights reserved. No part of this Manual may be displayed, reproduced or distributed in any form or by any means, without the prior written permission of the publisher, or used beyond the limited distribution to teachers and educators permitted by McGraw-Hill for their individual course preparation. If you are a student using this Manual, you are using it without permission.

4

f '( x )  f ( xi 1 )

2 x  ( xi  h)  ( xi  h) h h  f ( xi ) i  f ( xi 1 ) h (  h) 2h ( h) h(2h) f ( xi 1 ) f ( xi1 )  f ( xi1 )  f ( xi 1 )  0 

2h 2h 21.8 f uncti on [d, ea, i t er ] =r ombdi f f ( f unc, x, es, maxi t , varargi n) % r omber g: Romber g i nt egr at i on quadr at ure % q = r omber g( f unc, a, b, es, maxi t , var ar gi n) : % Romberg i nt egr ati on. % i nput : % f unc = name of f unct i on t o be i nt egr ated % a, b = i nt egr at i on l i mi t s % es = desi r ed r el at i ve er r or ( def aul t = 0. 000001%) % maxi t = maxi mum al l owabl e i t er at i ons ( def aul t = 50) % p1, p2, . . . = addi t i onal parameter s used by f unc % out put : % q = i nt egr al est i mate % ea = appr oxi mat e r el at i ve err or ( %) % i t er = number of i t er at i ons

2h

i f nar gi n<2, er r or ( ' at l east 2 i nput ar gument s requi r ed' ) , end i f nar gi n<3| i sempt y( es) , es=0. 000001; end i f nar gi n<4| i sempt y( maxi t ) , maxi t =50; end n = 1; DY( 1, 1) = dydxnew( f unc, x, n, var argi n{: }) ; i t er = 0; ea=100; whi l e i t er
Test the program by evaluating the deri vative of f ( x) = e –0.5 x at x = 1. The exact result is f '(1) = –0.5e –0.5(1) = –0.30326533. >> f or mat l ong >> f =@( x) exp( - 0. 5*x) ; >> [ d, ea, i t er ] =r o mbdi f f ( f , 1) d = - 0. 30326532985552 ea = 7. 603832401115816e- 008 i t er = 3


5

21.9 f unct i on dy=di f f uneq( x, y) % di f f uneq: di f f er ent i at i on of unequal l y- spaced dat a % dy=di f f uneq( f unc, x, es, maxi t , varargi n) : % i nput : % x = vect or of i ndependent var i abl e % y = vect or of dependent var i abl e % out put : % dy = vect or of der i vat i ve est i mates i f nar gi n<2, er r or ( ' at l eas t 2 i nput argument s r equi r ed' ) , end n=l engt h( x); f or i = 1: n dy( i ) = dyuneq( x, y, n, x( i ) ) ; end end f unct i on dydx=dyuneq( x, y, n, xx) i f xx <= x( 2) dydx = DyDx( xx, x( 1) , x( 2) , x( 3) , y( 1) , y( 2) , y( 3) ) ; el sei f xx >= x( n - 1) dydx = DyDx(xx, x(n - 2), x(n - 1), x( n) , y(n - 2) , y(n - 1) , y(n) ) ; el se f or i i = 2: n - 2 i f xx >= x( i i ) & xx <= x( i i + 1) i f xx - x( i i ) < x( i i +1) - xx %i f t he unknown i s cl oser t o t he l ower end of t he r ange, %x(i i ) wi l l be chosen as t he mi ddl e poi nt dydx = DyDx( x x, x( i i - 1) , x( i i ) , x( i i + 1) , y( i i - 1) , . . . y( i i ) , y( i i + 1) ) ; el s ei f xx - x( i i - 1) == x( i i +1) - xx & x( i i ) - x( i i - 1)
The program can be applied to the data with the following script: x=[ 0. 6 1. 5 1. 6 2. 5 3. 5] ; f x=[ 0. 9036 0. 3734 0. 3261 0. 08422 0. 01596] ; dydx=di f f uneq( x, f x) dydx = PROPRIETARY MATERIAL. © The McGraw-Hill Companies, Inc. All rights reserved. No part of this Manual may be displayed, reproduced or distributed in any form or by any means, without the prior written permission of the publisher, or used beyond the limited distribution to teachers and educators permitted by McGraw-Hill for their individual course preparation. If you are a student using this Manual, you are using it without permission.

6 - 0. 6936

- 0. 4846

- 0. 4526

- 0. 1738

0. 0373

The results can be compar ed with the true deriv atives which can be calculated with the analytical solution, x x f ( x) =5e –2 – 10 xe –2 . The results can be displayed graphically below where the computed values are represented as points and the true values as the curve. dyt r ue=5*exp( - 2*x)- 10*x. *exp( - 2*x) xpl ot =[ 0: 0. 125: 4] ; dyt r uepl ot =5*exp( - 2*xpl ot ) - 10*xpl ot . *exp( - 2*xpl ot ) ; pl ot ( x, dydx, ' o' , xpl ot , dyt r uepl ot ) l egend( ' numer i cal ' , ' anal yt i cal ' , ' l ocat i on' , ' bes t ' ) dyt r ue = - 0. 3012

- 0. 4979

- 0. 4484

- 0. 1348

- 0. 0274

6 numerical analytical

4 2 0 -2 0

1

2

3

4

As can be seen, the results leave something to be desired, particularly at the left end of the interval. The poor performance is in part due to the highly irregular spacing of the first three points. 21.10 f uncti on [ dydx, d2ydx2] = di f f eq( x, y) n = l engt h( x ) ; i f l engt h( y) ~=n, err or( ' x and y must be same l engt h' ) , end i f any( di f f ( di f f ( x) ) ~=0) , er r or ( ' unequal s pac i ng' ) , end i f l engt h( x) <4, er r or ( ' at l eas t 4 val ues r e qui r e d' ) , end dx=x(2)- x(1); f or i =1: n i f i ==1 dydx( i ) =( - y(i +2) +4*y( i +1) - 3*y( i ) ) / d x / 2 ; d2ydx2( i ) =( - y( i +3) +4*y( i +2) - 5*y( i +1) +2*y ( i ) ) / dx^2; el s ei f i ==n dydx( i ) =( 3* y( i ) - 4* y( i - 1) +y( i - 2) ) / dx/ 2; d2ydx2( i ) =( 2*y( i ) - 5*y( i - 1) +4*y( i - 2) - y( i - 3) ) / dx^2; el se dydx( i ) =( y( i +1) - y( i - 1) ) / dx/ 2; d2ydx2( i ) =( y( i +1) - 2*y( i ) +y( i - 1) ) / dx^2; end end subpl ot ( 2, 1, 1) ; pl ot ( x, dydx) ; gr i d; t i t l e( ' Fi r s t der i vat i ves ' ) subpl ot ( 2, 1, 2) ; pl ot ( x, d2ydx2) ; gr i d; t i t l e( ' Sec ond der i vat i ves ' )

The M-file can be run for the data from Prob. 21.11: >> t =[ 0 25 50 75 100 125] ; >> y=[ 0 32 58 78 92 100] ; >> [ dydx, d2ydx2] = di f f eq( t , y) dydx = PROPRIETARY MATERIAL. © The McGraw-Hill Companies, Inc. All rights reserved. No part of this Manual may be displayed, reproduced or distributed in any form or by any means, without the prior written permission of the publisher, or used beyond the limited distribution to teachers and educators permitted by McGraw-Hill for their individual course preparation. If you are a student using this Manual, you are using it without permission.

7 1. 4000 d2ydx2 = - 0. 0096

1. 1600

0. 9200

0. 6800

0. 4400

0. 2000

- 0. 0096

- 0. 0096

- 0. 0096

- 0. 0096

- 0. 0096

21.11 The first forward difference formula of O(h2) from Fig. 21.3 can be used to estimate the velocity for the first point at t = 0, f '(0) 

58  4(32)  3(0) 2(25)

 1.4

km s

The acceleration can be estimated with the second forward difference formula of O(h2) from Fig. 21.3 f "(0) 

78  4(58)  5(32)  2(0) (25)2

  0.0096

km s2

For the interior points, centered difference formulas of O(h2) from Fig. 21.5 can be used to estimate the velocities and accelerations. For example, at the second point at t = 25, f '(25)  f "(25) 

58  0 2(25)

 1.16

km

58  2(32)  0 (25)2

s

  0.0096

km s2

For the final point, backward difference formulas of O(h2) from Fig. 21.4 can be used to estimate the velocities and accelerations. The results for all values a re summarized in the following table. t 0 25

y 0 32

v 1.40 1.16

a -0.0096 -0.0096


8 50 75 100 125

58 78 92 100

0.92 0.68 0.44 0.20

-0.0096 -0.0096 -0.0096 -0.0096

21.12 Although the first three points are equally spaced, the remaining values are unequally spaced. Therefore, a good approach is to use Eq. 21.21 to perform the differentiation for all points. The results are summarized below: t 0 0.52 1.04 1.75 2.37 3.25 3.83

x 153 185 208 249 261 271 273

v 70.19231 52.88462 49.94473 37.25169 16.05181 6.59273 0.30382

a -47.09922 -19.46883 -10.82145 -26.58748 -24.50300 -10.80561 -10.88031

21.13 (a) x (t )  x( t i 1 ) 7.3  5.1 dx m v  x(t i )  i 1   0.55 dt 2h 4 s a 

d2x dt 2

 x(t i ) 

x (ti 1 )  2 x (ti )  x(t i 1 ) h2



7.3  2(6.3)  5.1 22

 0.05

m s2

(b) v a 

 x(ti  2 )  4 x(ti 1 )  3x (t i )



8  4(7.3)  3(6.3)

 0.575

m

2h 4 s  x(ti  3 )  4 x(ti  2 )  5x(ti 1 )  2 x(t i ) 8.4  4(8)  5(7.3)  2(6.3) h



2

2

2

 0.075

m s2

(c) v a 

3 x(ti )  4x (ti 1 )  x(t i 2 )



3(6.3)  4(5.1)  3.4

 0.475

m

2h 4 s 2 x (ti )  5 x (ti 1 )  4x (ti 2 )  x(t i 3 ) 2(6.3)  5(5.1)  4(3.4)  1.8 h



2

2

2

 0.275

m s2

21.14 d   (ti 1 )   ( t i 1 ) 0.67  0.70       0.0075 rad/s dt 2h 4 dr r (ti 1 )  r( ti  1 ) 6030  5560  r    117.5 m/s dt 2h 4   

 r 



v 

d 2 dt

2

d 2 r dt

2



 (ti 1 )  2 (ti )   (t i 1 ) h



r (ti 1 )  2r (ti )  r (ti  1 )

 117.5 er

a 

2

h

2



0.67  2(0.68)  0.70 (2)



2

6030  2(5800)  5560 (2)

2

 0.0025 rad/s2  2.5 m/s2

 43.5 e   2.82625 er 12.7375 e





21.15 The following script is developed to solve the problem and generate the plot: PROPRIETARY MATERIAL. © The McGraw-Hill Companies, Inc. All rights reserved. No part of this Manual may be displayed, reproduced or distributed in any form or by any means, without the prior written permission of the publisher, or used beyond the limited distribution to teachers and educators permitted by McGraw-Hill for their individual course preparation. If you are a student using this Manual, you are using it without permission.

9

cl ear , cl c, c l f t =[ 1 2 3. 25 4. 5 6 7 8 8. 5 9. 3 10] ; v=[ 10 12 11 14 17 16 12 14 14 10] ; f or i = 2: 4 p=pol yf i t ( t , v, i ) ; dvdt =pol yder( p) ; f or j = 1: l engt h( t ) vp=pol yval ( dvdt , t ) ; end s ubpl ot ( 3, 1, i - 1) , pl ot ( t , vp) , gr i d end 5 0 -5 0 5

2

4

6

8

10

2

4

6

8

10

2

4

6

8

10

0 -5 0 5 0 -5 0

21.16 >> x=- 2: . 1: 2; >> f =@( x) 1/ sqrt ( 2*pi ) *e xp( - ( x. ^2) / 2) ; >> y=f ( x) ; >> d=di f f ( y) . / di f f ( x ) ; >> x=- 1. 95: . 1: 1. 95; >> d2=di f f ( d) . / di f f ( x) ; >> x=- 1. 9: . 1: 1. 9; >> pl ot ( x, d2) ; gr i d


10

Thus, inflection points (d 2 y/dx2 = 0) occur at –1 and 1. 21.17 >> x=[ - 2 - 1. 5 - 1 - 0. 5 0 0. 5 1 1. 5 2] ; >> y=[ 0. 05399 0. 12952 0. 24197 0. 35207 0. 39894 0. 35207 0. 24197 0. 12952 0. 05399] ; >> d=di f f ( y) . / di f f ( x) ; >> x=- 1. 75: . 5: 1. 75; >> d2=di f f ( d) . / di f f ( x) ; >> x=- 1. 5: . 5: 1. 5; >> pl ot ( x, d2) ; gr i d

Thus, inflection points (d 2 y/dx2 = 0) occur at –1 and 1. 21.18 PROPRIETARY MATERIAL. © The McGraw-Hill Companies, Inc. All rights reserved. No part of this Manual may be displayed, reproduced or distributed in any form or by any means, without the prior written permission of the publisher, or used beyond the limited distribution to teachers and educators permitted by McGraw-Hill for their individual course preparation. If you are a student using this Manual, you are using it without permission.

11 %Cent ered Fi ni t e Di f f erence Fi r st & Second Deri vat i ves of Or der O( dx^2) %Usi ng di f f ( y) dx=1. ; y=[ 1. 4 2. 1 3. 3 4. 8 6. 8 6. 6 8. 6 7. 5 8. 9 10. 9 10] ; dyf =di f f ( y) ; % Fi r st Der i vat i ve Cent er ed FD usi ng di f f n=l engt h( y); f or i =1: n- 2 dydxc(i ) =( dyf ( i +1) +dyf ( i ) ) / ( 2*dx) ; end %Second Der i vat i ve Cent er ed FD usi ng di f f dy2dx2c=di f f ( dyf) / ( dx*dx); f pr i nt f ( ' f i r s t der i vat i ve \ n' ) ; f pr i nt f ( ' %f \ n' , dydxc) f pr i nt f ( ' s ec ond der i vat i ve \ n' ) ; f pr i nt f ( ' %f \ n' , dy2dx2c ) Results: f i r s t der i vat i ve 0. 950000 1. 350000 1. 750000 0. 900000 0. 900000 0. 450000 0. 150000 1. 700000 0. 550000 second deri vat i ve 0. 500000 0. 300000 0. 500000 - 2. 200000 2. 200000 - 3. 100000 2. 500000 0. 600000 - 2. 900000 21.19 Parts (a) through (d) can be solved with the following script: % % % % % %

Fi ni t e Di f f er ence Appr oxi mat i on of sl ope For f ( x) =exp( - 2*x) - x f ' ( x) =- 2*exp( - 2*x) - 1 Cent er ed di f f . df / dx=( f ( i +1) - f ( i - 1) ) / 2dx + O( dx^2) Fwd. di f f . df / dx=( - f ( i +2) +4f ( i +1) - 3f ( i ) ) / 2dx + O( dx^2) Bkwd. di f f . df / dx=( 3f ( i ) - 4f ( i - 1) +f ( i - 2) ) / 2dx + O( dx^2)

x=2; f x=exp( - 2*x)- x; df dx2=- 2*exp( - 2*x) - 1; %appr oxi mat i on dx=0. 5: - 0. 01: . 01; f or i =1: l engt h( dx) %x- val ues at i +- dx and +- 2dx xp( i ) =x+dx(i ) ; x2p( i ) =x+2*dx( i ) ; xn( i ) =x- dx( i ) ; x2n( i ) =x- 2*dx( i ) ; %f ( x) - val ues at i +- dx and +- 2dx f p( i ) =exp( - 2* xp( i ) ) - xp( i ) ; f 2p( i ) =exp( - 2* x2p( i ) ) - x2p( i ) ; f n( i ) =exp( - 2* xn( i ) ) - xn( i ) ; PROPRIETARY MATERIAL. © The McGraw-Hill Companies, Inc. All rights reserved. No part of this Manual may be displayed, reproduced or distributed in any form or by any means, without the prior written permission of the publisher, or used beyond the limited distribution to teachers and educators permitted by McGraw-Hill for their individual course preparation. If you are a student using this Manual, you are using it without permission.

12 f 2n( i ) =exp( - 2* x2n( i ) ) - x2n( i ) ; %Fi ni t e Di f f . Appr oxi mat i ons Cdf dx( i ) =( f p( i ) - f n( i ) ) / ( 2* dx( i ) ) ; Fdf dx( i ) =( - f 2p( i ) +4* f p( i ) - 3* f x ) / ( 2* dx( i ) ) ; Bdf dx( i ) =( 3* f x- 4* f n ( i ) +f 2n( i ) ) / ( 2* dx( i ) ) ; end dx0=0; pl ot ( dx, Fdf dx, ' - - ' , dx, Bdf dx, ' - . ' , dx, Cdf dx, ' - ' , dx0, df dx2, ' *' ) gr i d t i t l e( ' For ward, Backward, and Cent ered Fi ni t e Di f f erence appr oxi mat i on Or der Corr ect ' ) xl abel ( ' Del t a x' ) yl abel ( ' df / dx' ) gt ext ( ' Cent er ed' ) ; gt ext ( ' For war d' ) ; gt ext ( ' Backwar d' )

2nd

21.20 The flow rate is equal to the derivative of volume with respect to time. Equation (21.21) can be used to compute the derivative as x0 = 1 x1 = 5 x2 = 8 f '(7)  1

f ( x0) = 1 f ( x1) = 8 f ( x2) = 16.4

2(7)  5  8 (1  5)(1  8)

8

2(7)  1  8 (5  1)(5  8)

 16.4

2(7)  1  5 (8  1)(8  5)

 0.035714  3.33333  6.247619  2.95

Therefore, the flow is equal to 2.95 cm3/s. 21.21 The velocity at the surface can be computed with Eq. (21.21) as x0 = 0 x1 = 0.002 x2 = 0.006

f ( x0) = 0 f ( x1) = 0.287 f ( x2) = 0.899


13

f '(0)  0

2(0)  0.002  0.006

 0.287

(0  0.002)(0  0.006)

2(0)  0  0.006 (0.002  0)(0.002  0.006)

 0.899

2(0)  0  0.002 (0.006  0)(0.006  0.002)

 0  215.25  74.9167  140.3333 Therefore, the shear stress can be computed as   1.8  10

5 N  s

1 N 140.3333  0.00253 2 s m m 2

21.22 Equation (21.21) can be used to compute the derivative as x0 = 0 x1 = 1 x2 = 3

f ( x0) = 0.06 f ( x1) = 0.32 f ( x2) = 0.6

f '(0)  0.06

2(0)  1  3 (0  1)(0  3)

 0.32

2(0)  0  3 (1  0)(1  3)

 0.6

2(0)  0  1 (3  0)(3  1)

 0.08  0.48  0.1  0.3

The mass flux can be computed as Mass flux  1.52  106

cm2 s

 0.3  106

g cm4

 4.56  1013

g cm2s

where the negative sign connotes transport from the sediments into the lake. The amount of mass transported into the lake can be computed as Mass transport  4.56 10

g

13

6

2

cm s

2

(3.6 10 m )

365 d 86,400 s 10,000 cm yr

d

m

2

2

kg 1000 g

 517.695 kg

21.23 For the equally-spaced points, we can use the second-order formulas from Figs. 21.3 through 21.5. For example, for the first point ( t = 0), we can use f '(0) 

0.77  4(0.7)  3(0.4) 2(10)

 0.0415

barrels min

However, the points around t = 30 are unequally spaced so we must use Eq. (21.21) to compute the derivative as x0 = 20 x1 = 30 x2 = 45 f '(30)  0.77

f ( x0) = 0.77 f ( x1) = 0.88 f ( x2) = 1.05

2(30)  30  45 (20  30)(20  45)

 0.88

2(30)  20  45 (30  20)(30  45)

 1.05

2(30)  20  30 (45  20)(45  30)

 0.011133

All the results can be summarized as t 0 10 20 30

V 0.4 0.7 0.77 0.88

Derivative 0.0415 0.0185 0.009 0.011133

Method Forward Centered Centered Unequal


14 45 60 75

1.05 1.17 1.35

0.009667 0.01 0.014

Centered Centered Backward

The derivatives can be plotted versus time as 1.6

0.05 0.04

1.2

0.03 0.8 0.02 0.4

0.01

0

0 0

10

20

30

40 V

50

60

70

80

Derivative

21.24 First printing: Because the values are equally spaced and the changes are equal, we can use a first-order approximation to determine the slope at x = 0 as x0 = 0 x1 = 0.08 f '(0) 

f ( x0) = 19 f ( x1) = 17

17  19 0.08

  25

C m

The coefficient of thermal conductivity can then be estimated as k = –60 W/m2/(–25oC/m) = 2.4 W/( oCm). Second (and later) printing: To make the problem more interesting, the surface value is changed to 20.2. For this case, the second-order forward difference from Fig. 21.3 is used to compute the derivative as x0 = 0 x1 = 0.08 x2 = 0.16 f '(0) 

f ( x0) = 20.2 f ( x1) = 17 f ( x2) = 15

15  4(17)  3(20.2) 0.16

  47.5

C m

The coefficient of thermal conductivity can then be estimated as k = –60 W/m2/(–47.5oC/m) = 1.5 W/(oCm). 21.25 First, we can estimate the areas by numerically differentiating the volume data. Because the values are equally spaced, we can use the second-order difference formulas from Fig. 21.3 to compute the derivatives at each depth. For example, at the first depth, we can use the forward difference to compute As (0)  

dV dz

(0)  

1,963,500  4(5,105,100)  3(9,817,500) 8

 1,374, 450 m2

For the interior points, second-order centered differences can be used. For example, at the second point at ( z = 4 m), PROPRIETARY MATERIAL. © The McGraw-Hill Companies, Inc. All rights reserved. No part of this Manual may be displayed, reproduced or distributed in any form or by any means, without the prior written permission of the publisher, or used beyond the limited distribution to teachers and educators permitted by McGraw-Hill for their individual course preparation. If you are a student using this Manual, you are using it without permission.

15

As (4)  

dV

(4)  

dz

1,963,500  9,817,500 8

 981, 750 m2

The other interior points can be determined in a similar fashion As (8)  

dV dz

As (12)  

(8)  

dV dz

392,700  5,105,100 8

(12)  

0  1, 963, 500 8

 589, 050 m 2

 245, 437.5 m 2

For the last point, the second-order backward formula yields As (16)  

dV dz

(16)  

3(0)  4(392,700)  1,963,500 8

 49, 087.5 m2

Since this is clearly a physically unrealistic result, we will assume that the bottom area is 0. The results are summarized in the following table along with the other quantities needed to determine the average concentration. 3

z, m 0 4 8 12 16

V, m 9817500 5105100 1963500 392700 0

c, g/m 10.2 8.5 7.4 5.2 4.1

3

2

A s , m 1374450.0 981750.0 589050.0 245437.5 0

c A s 14019390 8344875 4358970 1276275 0

The necessary integrals can then be evaluated with the multi-segment Simpson’s 1/3 rule,

 

z

As ( z ) dz  (16  0)

0 z

1, 374, 450  4(981, 750  245, 437.5)  2(589, 050)  0

c ( z ) As ( z ) dz  (16  0)

0

 9, 948, 400 m3

12 14, 0 19, 390  4(8, 344, 875  1, 2 76, 2 75)  2( 4, 358, 970)  0 12

 81,629,240 g

The average concentration can then be computed as

c 



Z

c( z ) As ( z ) dz

0



Z

0

As ( z ) dz



81,629,240 9,948,400

 8.205263

g m

3

21.26 For the equispaced data, we can use the second-order finite divided difference formulas from Figs. 21.3 through 21.5. For example, for the first point, we can use di dt



0.32  4(0.16)  3(0)  1.6 0.2  0

For the intermediate equispaced points, we can use centered differences. For example, for the second point di dt



0.32  0 0.2  0

 1.6


16

For the last point, we can use a backward difference di dt



3(2)  4(0.84)  0.56 0.7  0.3

8

One of the points ( t = 0.3) has unequally spaced neighbors. For this point, we can use Eq. (23.9) to compute the derivative as x0 = 0.2 f ( x0) = 0.32 x1 = 0.3 f ( x1) = 0.56 x2 = 0.5 f ( x2) = 0.84 f '(0)  0.32

2(0.3)  0.3  0.5 (0.2  0.3)(0.2  0.5)

 0.56

2(0.3)  0.2  0.5 (0.3  0.2)(0.3  0.5)

 0.84

2(0.3)  0.2  0.3 (0.5  0.2)(0.5  0.3)

 2.066667

We can then multiply these derivative estimates by the inductance. All the results are summarized below: t 0 0.1 0.2 0.3 0.5 0.7

i 0 0.16 0.32 0.56 0.84 2

di /dt 1.6 1.6 2 2.066667 3.6 8

VL 6.4 6.4 8 8.266667 14.4 32

21.27 We can solve Faraday’s law for inductance as



L 

t

0

V L dt i

We can evaluate the integral using a combination of the trapezoidal and Simpson’s rules, I  (20  0)

0  4(18)  29 6

 (180  120)

 (80  20)

35  26

2

29  3(44  49)  46

 (280  180)

8 26  15

 (120  80)

46  35

2 15  7

 (400  280) 2 2  336.6667  2655  1620  1830  2050  1320  9811.667

The inductance can then be computed as L 

9811.667 2

 4905.833

volts ms

s

A

1000 ms

 4.905833 H

21.28 Because the data is equispaced, we can use the second-order finite divided difference formulas from Figs. 21.3 through 21.5. For the first point, we can use dT dt



30  4(44.5)  3(80)  9.2 10  0

For the intermediate points, we can use centered differences. For example, for the second point


17

dT

30  80



dt

10  0

 5

We can analyze the other interior points in a similar fashion. For the last point, we can use a backward difference dT

3(20.7)  4(21.7)  24.1



dt

25  15

 0.06

All the values can be tabulated as t 0 5 10 15 20 25

T 80 44.5 30 24.1 21.7 20.7

dT/dt -9.2 -5 -2.04 -0.83 -0.34 -0.06

T Ta 60 24.5 10 4.1 1.7 0.7

If Newton’s law of cooling holds, we can plot dT /dt versus T – T a and the points can be fit with a linear regression with zero intercept to estimate the cooling rate. As in the following plot, the result is k = 0.1618/min. 0

10

20

30

40

50

60

0 -2 -4 -6 -8

y = -0.1618x

-10

R = 0.9758

2

21.29 The following script solves the problem. Note that the derivative is calculated with a centered difference, dV dT



V450 K  V 350 K

100 K

The following script evaluates the derivative with centered finite differences and the integral with the t r apz function, V=[ 220 250 282. 5; 4. 1 4. 7 5. 23; 2. 2 2. 5 2. 7; 1. 35 1. 49 1. 55; 1. 1 1. 2 1. 24; . 9 . 99 1. 03; . 68 . 75 . 78; . 61 . 675 . 7; . 54 . 6 . 62] ; P=[ 0. 1 5 10 20 25 30 40 45 50] ' ; T=[ 350 400 450] ' ; n=l engt h( V) ; dVdt =( V( 1: n, 3) - V( 1: n, 1) ) / ( T( 3) - T( 1) ) ; i nt egr and=V( 1: n, 2) - T( 2) *dVdt ; H=t r apz( P, i nt egr and)

When the script is run the result is H = 21.4410. The following table displays all the results


18

P, atm 0.1 5 10 20 25 30 40 45 50

T = 350K 220 4.1 2.2 1.35 1.1 0.9 0.68 0.61 0.54

T = 400K 250 4.7 2.5 1.49 1.2 0.99 0.75 0.675 0.6

T = 450K 282.5 5.23 2.7 1.55 1.24 1.03 0.78 0.7 0.62

dV/dT 0.625 0.0113 0.005 0.002 0.0014 0.0013 0.001 0.0009 0.0008

T (dV/dT)p 0 0.18 0.5 0.69 0.64 0.47 0.35 0.315 0.28 Total Integral =

Trap

V

0.441 1.7 5.95 3.325 2.775 4.1 1.6625 1.4875 21.441

21.30 (a) First, the distance can be converted to meters. Then, Eq. (21.21) can be used to compute the derivative at the surface as x0 = 0 x1 = 0.01 x2 = 0.03

f ( x0) = 900 f ( x1) = 480 f ( x2) = 270

f '(0)  900

2(0)  0.01  0.03 (0  0.01)(0  0.03)

 480

2(0)  0  0.03 (0.01  0)(0.01  0.03)

 270

2(0)  0  0.01 (0.03  0)(0.03  0.01)

 52,500

K m

The heat flux can be computed as Heat flux  0.028

K  W  52,500   1, 470 2  smK  m m J

(b) The heat transfer can be computed by multiplying the flux by the area

Heat transfer  1, 470

W m2

(200 cm  50 cm)

m2 10,000 cm2

 1, 470 W

21.31 (a) The pressure drop can be determined by integrating the pressure gradient p 



x2



x1

8 Q  r ( x) 4

d x

After converting the units to meters, a table can be set up holding the data and the integrand. The trapezoidal and Simpson’s rules can then be used to integrate this data as shown in the last column of the table. x, m 0 0.02 0.04 0.05 0.06 0.07 0.1

r , m 0.002 0.00135 0.00134 0.0016 0.00158 0.00142 0.002

integrand -7957.75 -38333.2 -39490.3 -19428.1 -20430.6 -31315.3 -7957.75 Sum

integral

method

-1338.5392

Simp 1/3

-713.9319 -589.0959 -2641.5670

Simp 3/8 Trap

Therefore, the pressure drop is computed as –2,641.567 N/m 2.


19

(b) The average radius can also be computed by integration as

 r 

x2

r (x ) dx

x1

x2  x1

The numerical evaluations can again be determined by a combination of the trapezoidal and Simpson’s rules. x, m 0 0.02 0.04 0.05 0.06 0.07 0.1

r , m 0.00200 0.00135 0.00134 0.00160 0.00158 0.00142 0.00200 Sum

integral

method

5.83E-05

Simp 1/3

4.61E-05 5.13E-05 0.0001557

Simp 3/8 Trap

Therefore, the average radius is 0.0001557/0.1 = 0.001557 m. This value can be used to compute a pressure drop of

 p 

dp dx

x  

8  Q

 r 4

x  

8(0.005)0.00001

 (0.001557) 4

0.1  2166.95

N m

2

Thus, there is less pressure drop if the radius is at the constant mean value. (c) The average Reynolds number can be computed by first determining the average velocity as v

Q Ac

0.00001



 (0.001557)

2



0.00001 7.6152 10

6

 1.31317

m s

Then the Reynolds number can be computed as Re 

3

1 10 (1.31317)0.0031138 0.005

 817.796

21.32 Using an approach based on Example 21.4, the following script determines first through fourth derivative estimates by taking differences of differences. A plot of the results is also generated. T=[ 300 400 500 600 700 800 900 1000] ; Cp=[ 82. 888 112. 136 136. 933 157. 744 175. 036 189. 273 200. 923 210. 450] ; dCp=di f f ( Cp) . / di f f ( T) ; n=l engt h( T); Tm=( T( 1: n- 1) +T( 2: n) ) . / 2; dCp2=di f f ( dCp) . / di f f ( Tm) ; n2=l engt h( Tm) ; Tm2=( Tm( 1: n2- 1) +Tm( 2: n2) ) . / 2; dCp3=di f f ( dCp2) . / di f f ( Tm2) ; n3=l engt h( Tm2) ; Tm3=( Tm2( 1: n3- 1) +Tm2( 2: n3) ) . / 2; dCp4=di f f ( dCp3) . / di f f ( Tm3) ; n4=l engt h( Tm3) ; Tm4=( Tm3( 1: n4- 1) +Tm3( 2: n4) ) . / 2; PROPRIETARY MATERIAL. © The McGraw-Hill Companies, Inc. All rights reserved. No part of this Manual may be displayed, reproduced or distributed in any form or by any means, without the prior written permission of the publisher, or used beyond the limited distribution to teachers and educators permitted by McGraw-Hill for their individual course preparation. If you are a student using this Manual, you are using it without permission.

20 subpl ot ( 4, 1, 1) , pl ot ( Tm, dCp) gr i d, xl i m( [ mi n( T) max( T) ] ) , t i t l subpl ot ( 4, 1, 2) , pl ot ( Tm2, dCp2) gr i d, xl i m( [ mi n( T) max( T) ] ) , t i t l subpl ot ( 4, 1, 3) , pl ot ( Tm3, dCp3) gr i d, xl i m( [ mi n( T) max( T) ] ) , t i t l subpl ot ( 4, 1, 4) , pl ot ( Tm4, dCp4) gr i d, xl i m( [ mi n( T) max( T) ] ) , t i t l

e( ' f i r s t der i vat i ve' ) e( ' s ec ond der i vat i ve' ) e( ' t h i r d der i vat i ve' ) e( ' f our t h der i vat i ve' )

The fact that (1) the second derivative appears linear, (2) the third derivative is close to being constant, and (3) the fourth derivative oscillates in a random fashion close to zero suggests that the data was generated with a fourth-order polynomial (which in fact was the case). 21.33 Equation (21.21) can be used to estimate the derivatives at each temperature. This can be done using the M-file developed in Prob. 21.9 with the result >> T=[ 750 800 900 1000] ; >> h=[ 29629 32179 37405 42769] ; >> cp=di f f uneq( T, h) cp = 50. 5800

51. 4200

52. 9500

54. 3300

These results, which have units of kJ/(kmol K), can b e converted to the desired units as in 3 3  J  mol 10 g  kJ  kmol 10 J c p       cp    kmol K  103 mol kJ 12.011  2(15.9994)  g kg  kg K 

>> cp=cp* 1e3/ ( 12. 011+2*15. 9994) cp = 1. 0e+003 * 1. 1493

1. 1684

1. 2031

1. 2345


21

21.34 The following script determines dc/dt using an M-file that implements Eq. (21.21). Linear regression is then used to fit the log transformed equation. The best-fit parameters are then used to compute the model coefficients n and k with the following script: t =[ 0 5 15 30 45] ; c=[ 0. 750 0. 594 0. 420 0. 291 0. 223] ; dcdt =di f f uneq( t , c) p=pol yf i t ( l og10( c ) , l og10( - dc dt ) , 1) n=p( 1) k=10^p( 2)

>> pr ob1934 dcdt = - 0. 0358 p = 2. 1419 n = 2. 1419 k = 0. 0773

- 0. 0266

- 0. 0139

- 0. 0066

- 0. 0025

- 1. 1121

21.35 Finite-divided differences can be used to numerically estimate the derivatives. For the first time ( t = 0), a forward difference of O(h2) can be employed to give do dt

(0) 

3(10)  4(7.11)  4.59 0.25

  24.6

mg L d

Centered differences of O(h2) can be used for the interior data. For example, for the second interval, do dt

(0.125) 

4.59 10 0.25

  21.64

mg Ld

A backward difference of O(h2) can be employed for the last time ( t = 0.75) to give do dt

(0.75) 

3(0.03)  4(0.33)  1.15 0.25

  0.32

mg Ld

The SODs can then be computed with the formula from the problem statement. For example, for the first point, do g SOD   H  0.1(24.6)  2.46 dt m2 d

All the results can be tabulated along with the SOD: t (d) 0 0.125 0.25 0.375 0.5 0.625 0.75

o (mg/L) 10.00 7.11 4.59 2.57 1.15 0.33 0.03

do /dt -24.60 -21.64 -18.16 -13.76 -8.96 -4.48 -0.32

SOD 2.460 2.164 1.816 1.376 0.896 0.448 0.032


22

The plots can be developed as 3

SOD (g/m 2/d)

2

1

0 0

0.2

0.4

0.6

0.8

t (d) SOD (g/m 2/d)

3

2

1

0 0

2

4

6

8

10

o (mg/L)

21.36 (a) Displacement can be determined by integrating the slope equation (P21.35), y ( x) 

w0

120 EIL



x

0

5 x 4  6 L2 x2  L4 dx

Because this is a polynomial, the exact solution can be evaluated analytically as y 

w0

  x5  2 L2 x3  L4 x   0.0000115741x 5  0.000208333x3  0.0009375x  120 EIL 

Numerical solutions can be obtained with the trapezoidal. For example, the displacement at x = 0.125 can be computed as  (0) 

250000

 5(0) 4  2(3)2 (0)2  (3)4   9.375  104  120(2 10 )0.0003(3)  11

 (0.125) 

250000

 5(0.125)4  2(3)2 (0.125)2  (3)4   9.27749  104  120(2 10 )0.0003(3)  11

y (0.125)  (0.125  0)

9.375 104  9.27749 10 4 2

 0.000116578 m

At x = 0.25:  (0.25) 

250000

 5(0.25)4  2(3)2 (0.25)2  (3)4   8.98664  104  120(2  10 )0.0003(3)  11


23 9.27749 10   8.98664 10  4

y (0.25)  0.000116578  (0.25  0.125)

4

2

 0.000116578  0.000114151=  0.000230729

The remainder of the calculation along with the analytical solution is summarized in the following table and graph: x

y-analytical

0 0.125 0.25 0.375 0.5 0.625 0.75 0.875 1 1.125 1.25 1.375 1.5 1.625 1.75 1.875 2 2.125 2.25 2.375 2.5 2.625 2.75 2.875 3

-0.0009375 -0.0009277 -0.0008987 -0.0008508 -0.0007849 -0.0007022 -0.0006042 -0.0004929 -0.0003704 -0.0002392 -0.0001022 0.0000373 0.0001758 0.0003094 0.0004338 0.0005445 0.0006366 0.0007047 0.0007434 0.0007466 0.0007082 0.0006214 0.0004794 0.0002748 0.0000000

Trap Rule 0.0000000 -0.0001166 -0.0001142 -0.0001093 -0.0001022 -0.0000929 -0.0000817 -0.0000686 -0.0000540 -0.0000381 -0.0000213 -0.0000041 0.0000133 0.0000303 0.0000464 0.0000611 0.0000738 0.0000838 0.0000905 0.0000931 0.0000909 0.0000831 0.0000688 0.0000471 0.0000172

0.0000000 -0.0001168 -0.0002311 -0.0003407 -0.0004431 -0.0005362 -0.0006180 -0.0006867 -0.0007407 -0.0007789 -0.0008003 -0.0008044 -0.0007910 -0.0007606 -0.0007141 -0.0006527 -0.0005787 -0.0004946 -0.0004037 -0.0003102 -0.0002188 -0.0001352 -0.0000658 -0.0000179 0.0000000

y = Trap Rule 0.0000000 -0.0001166 -0.0002307 -0.0003401 -0.0004423 -0.0005352 -0.0006169 -0.0006855 -0.0007394 -0.0007775 -0.0007988 -0.0008029 -0.0007896 -0.0007593 -0.0007128 -0.0006517 -0.0005779 -0.0004940 -0.0004035 -0.0003104 -0.0002195 -0.0001364 -0.0000676 -0.0000204 -0.0000033

The following plot indicates good agreement between the analytical and numerical results. Note that we could improve the results by using a finer step size (i.e., a smaller value of  x), a more refined numerical integration formula (Simpson’s 1/3 rule) or a higher accuracy method (e.g., Romberg integration or the MATLAB quad function). Deflection vers us x 0.0000 -0.0002

0

1

2

3

-0.0004 -0.0006 Ana ly tic al Numerical

-0.0008 -0.0010

(b) The moment can also be evaluated analytically: M ( x)  EI

d  dx



w0

120 L

( 20 x3  12 L2 x)

Numerical solutions can be obtained with finite differences. For example, the moment at x = 0 can be computed with a forward finite difference estimate for the first derivative (Table 21.3).


24

M (0)  EI

d  dx

 2 1011 (0.0003)

3(9.375  104 )  4(9.277  104 )  (8.987  104 ) 2(0.125)

 40.69

Centered differences (Table 21.5) can then be used for the interior points and a backward difference (Table 21.4) for x = 3. For example, for x = 0.125, M (0.125)  2 1011 (0.0003)

8.987  104  (9.375  104 ) 2(0.125)

 9320.75

The remainder of the calculation along with the analytical solution is summarized in the following graph and the table at the end of this problem solution. Moment versus x

100000

0 0

0.5

1

-100000

1.5

2

2.5

3

Analy tic al Numerical

-200000

(c) The shear can also be evaluated analytically: V ( x)  EI

2

d  dx

2



w0

120 L

( 60 x 2  12 L2 )

Numerical solutions can be obtained with finite differences. For example, the shear at x = 0 can be computed with a forward finite difference estimate for the second derivative (Table 21.3) V (0)  2  1011 (0.0003)

(8.508  104 )  4(8.987  104 )  5(9.277  104 )  2(9.375  104 ) (0.125)2

 76193.58

Centered differences (Table 21.5) can then be used for the interior points and a backward difference (Table 21.4) for x = 3. For example, for x = 0.125, 11

V (0)  2  10 (0.0003)

8.987  104  2(9.277  104 )  9.375  104 (0.125)2

 74240.45

The remainder of the calculation along with the analytical solution is summarized in the following graph and the table below.


25

Shear versus x

100000 0 -100000

0

1

2

-200000

3

Analy tic al Numerical

-300000 -400000 x

M-anal M-num 0 0.0 40.7 0.125 9347.9 9320.7 0.25 18533.0 18478.7 0.375 27392.6 27311.2 0.5 35763.9 35655.4 0.625 43484.2 43348.5 0.75 50390.6 50227.9 0.875 56320.5 56130.6 1 61111.1 60894.1 1.125 64599.6 64355.5 1.25 66623.3 66352.0 1.375 67019.3 66720.9 1.5 65625.0 65299.5 1.625 62277.6 61924.9 1.75 56814.2 56434.5 1.875 49072.3 48665.4 2 38888.9 38454.9 2.125 26101.3 25640.2 2.25 10546.9 1 0058.6 2.375 -7937.3 -8452.7 2.5 -29513.9 -30056.4 2.625 -54345.7 -54915.4 2.75 -82595.5 -83192.3 2.875 -114426.0 -115049.9 3 -150000.0 -148738.6

V-anal 75000.0 74349.0 72395.8 69140.6 64583.3 58724.0 51562.5 43099.0 33333.3 22265.6 9895.8 -3776.0 -18750.0 -35026.0 -52604.2 -71484.4 -91666.7 -113151.0 -135937.5 -160026.0 -185416.7 -212109.4 -240104.2 -269401.0 -300000.0

V-num 76193.6 74240.5 72287.3 69032.1 64474.8 58615.5 51454.0 42990.5 33224.8 22157.1 9787.3 -3884.5 -18858.5 -35134.5 -52712.7 -71592.9 -91775.2 -113259.5 -136046.0 -160134.5 -185525.2 -212217.9 -240212.7 -269509.5 -298806.4

The above results indicate good agreement between the analytical and numerical results. We could improve the results by using a finer step size (i.e., a smaller value of  x). 21.37 From Prob. 21.36, dy dx

  ( x)

d dx



M ( x)

dM

EI

dx

 V ( x)

dV dx

  w( x)

The following script calls the function to generate and p lot the derivatives. Script: cl ear , cl c, cl f x=[ 0 0. 375 0. 75 1. 125 1. 5 1. 875 2. 25 2. 625 3] ; y=[ 0 - 0. 2571 - 0. 9484 - 1. 9689 - 3. 2262 - 4. 6414 - 6. 1503 - 7. 7051 - 9. 275] ; y=y/ 100; E=200e9; I =0. 0003; [ dydx, d2ydx2, d3ydx3, d4ydx4] =Deri vs( x, y) ; M=E*I *d2ydx2; V=E*I *d3ydx3; wx=E*I *d4ydx4; PROPRIETARY MATERIAL. © The McGraw-Hill Companies, Inc. All rights reserved. No part of this Manual may be displayed, reproduced or distributed in any form or by any means, without the prior written permission of the publisher, or used beyond the limited distribution to teachers and educators permitted by McGraw-Hill for their individual course preparation. If you are a student using this Manual, you are using it without permission.

26 subpl ot ( 5, 1, 1) pl ot ( x, y, ' o- ' ) subpl ot ( 5, 1, 2) pl ot ( x, dydx) subpl ot ( 5, 1, 3) pl ot ( x, M) subpl ot ( 5, 1, 4) pl ot ( x, V) subpl ot ( 5, 1, 5) pl ot ( x, wx) z=[ x' y' dydx M V wx] ' ; f pr i nt f ( ' x y dy/ dx M( x) V( x) f pr i nt f ( ' %10. 4f %10. 4f %10. 4f %10. 1f %10. 1f %10. 1f \ n' , z) ;

d4y/ dx4\ n' ) ;

The function, Der i vs , is used to determine the derivatives. Because the data are equi-spaced, we can use the formulas from Tables 21.3 through 21.5 to determine the derivatives with the following functions. Notice that the function Di f f 1 determines first derivatives and Di f f 2 determines second derivatives. f unct i on [ dydx, d2ydx2, d3ydx3, d4ydx4] =Deri vs( x, y) clf i f any( di f f ( di f f ( x) ) ~=0) , er r or ( ' unequal s pac i ng' ) , end i f any( di f f ( x) <=0) , er r or( ' not i n ascendi ng or der' ) , end dydx=Di f f Num( x, y) ; d2ydx2=Di f f Num2( x, y) ; d3ydx3=Di f f Num( x, d2ydx2) ; d4ydx4=Di f f Num( x, d3ydx3) ; f unct i on dydx=Di f f Num( x, y) n=l engt h( x); dydx=zeros( n) ; dx=x( 2) - x( 1) ; dydx(1)=( - 3*y(1)+4*y(2)- y( 3) ) / ( 2*dx) ; f or i = 2: n- 1 dydx( i ) =( y( i +1) - y( i - 1) ) / ( 2* dx) ; end dydx(n)=( 3*y(n) - 4*y(n- 1) +y(n- 2) ) / ( 2*dx); dydx=dydx( : , 1) ; end f unct i on d2ydx2=Di f f Num2( x, y) n=l engt h( x); d2ydx2=zer os( n) ; dx=x( 2) - x( 1) ; d2ydx2( 1) =( - y( 4) +4*y( 3) - 5*y( 2) +2*y( 1) ) / dx^2; f or i = 2: n- 1 d2ydx2( i ) =( y(i +1) - 2*y(i ) +y(i - 1) ) / dx^2; end d2ydx2( n) =( 2*y( n) - 5*y( n- 1) +4*y( n- 2) - y( n- 3) ) / dx^2; d2ydx2=d2ydx2( : , 1) ; end Results: x 0. 0000 0. 3750 0. 7500 1. 1250 1. 5000 1. 8750 2. 2500 2. 6250 3. 0000

y 0. 0000 - 0. 0026 - 0. 0095 - 0. 0197 - 0. 0323 - 0. 0464 - 0. 0615 - 0. 0771 - 0. 0927

dy/ dx - 0. 0011 - 0. 0126 - 0. 0228 - 0. 0304 - 0. 0356 - 0. 0390 - 0. 0408 - 0. 0417 - 0. 0421

M( x) - 2300586. 7 - 1852586. 7 - 1404586. 7 - 1010346. 7 - 673706. 7 - 399786. 7 - 195840. 0 - 64426. 7 66986. 7

V( x) 1194666. 7 1194666. 7 1122986. 7 974506. 7 814080. 0 637155. 6 447146. 7 350435. 6 350435. 6

d4y/ dx4 95573. 3 - 95573. 3 - 293546. 7 - 411875. 6 - 449801. 5 - 489244. 4 - 382293. 3 - 128948. 1 128948. 1


27

0 -0.05 -0.1

0

0.5

1

1.5

2

2.5

3

0.5

1

1.5

2

2.5

3

0.5

1

1.5

2

2.5

3

0.5

1

1.5

2

2.5

3

0.5

1

1.5

2

2.5

3

0 -0.02 -0.04 0 0 -1 -2 10 5 0 1 0 -1

x 10

6

0 x 10

5

0 x 10

0

6

21.38 (a)  f  3 y  3  3x 2  3(1)  3  3(1)2  3  x

 f  3 x  9 y 2  3(1)  9(1)2  6  y  f 3  xy (b)  f

 x

 



f ( x  x, y )  f ( x  x, y )

2x



f (1.0001,1)  f (0.9999,1)

2x 3

3

[3(1.0001)1  3(1.0001)  (1.0001)  3(1) ]  [3(0.9999)1 3(0.9999)  (0.9999)3  3(1)3 ] 0.0002 2.00029997  1.99969997 0.0002

 2.99999999

f ( x, y  y )  f ( x, y  y) f (1,1.0001)  f (1, 0.9999)  f   2y 2y  y





[3(1)1.0001  3(1)  (1)3  3(1.0001)3 ]  [3(1)0.9999  3(1)  (1)3  3(0.9999)3 ] 0.0002 1.99939991 2.00059991 0.0002

  6.00000003

f ( x  x, y  y )  f ( x  x, y  y )  f ( x  x, y  y )  f ( x  x, y  y)  2 f   xy 4xy



f (1.0001,1.0001)  f (1.0001,0.9999)  f (0.9999,1.0001)  f (0.9999,0.9999)

4 xy


28



1.99969991  2.00089985  1.99909985  2.00029991 4(0.0001)(0.0001)

 2.999999982

21.39 Script: cl ear , cl c, cl f f =@( x, y) exp( - ( x. ^2+y. ^2) ) ; [ x, y] =meshgr i d( - 3: 0. 25: 3, - 3: 0. 25: 3) ; z =f ( x, y) ; [ f x, f y] =gr a di ent ( z , 0. 25) ; cs=cont our ( x, y, z) ; cl abel ( cs) hol d on qui ver ( x, y, - f x, - f y) ; axi s s quar e hol d of f pause f =@( x, y) x. *exp( - ( x. ^2+y. ^2) ) ; [ x, y] =meshgr i d( - 3: 0. 25: 3, - 3: 0. 25: 3) ; z =f ( x, y) ; [ f x, f y] =gr a di ent ( z , 0. 25) ; cs=cont our ( x, y, z) ; cl abel ( cs); hol d on hol d on qui ver ( x, y, - f x, - f y) ; axi s s quar e hol d of f Output:

3

2 0.1 0.3 0.5

1

0.7 0.9 1

0

0.8 0.6 0.4

-1

0.2

-2

-3 -3

-2

-1

0

1

2

3


29

3

2

1 0.1 0.3

0

0.4

-0.4 -0.2 0 -0.3

0.2 -0.1

-1

-2

-3 -3

-2

-1

0

1

2

3

21.40 Script: cl ear , cl c, cl f [ x, y] =meshgr i d( - 3: 0. 25: 3, - 3: 0. 25: 3) ; z=peaks( x, y); [ f x, f y] =gr a di ent ( z , 0. 25) ; cs=cont our ( x, y, z) ; cl abel ( cs) hol d on qui ver ( x, y, - f x, - f y) ; axi s s quar e hol d of f Result:

3

2

4 6

1 2 -2

0

0

2

-1 -6

0

-4

-2

-2 0

-3 -3

-2

-1

0

1

2

3


sm ch (21).pdf

Recommend Documents