1 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. CHAPTER 28 28.1 The solution with the 2nd-order RK (Heun without corrector) can be laid out as t c k1 c k2 0 10 2 30 1 1.510 25 1.25 37.5 0.625 0.937520 34.375 0.78125 42.1875 0.390625 0.58593830 40.23438 0.488281 45.11719 0.244141 0.36621140 43.89648 0.305176 46.94824 0.152588 0.22888250 46.1853 0.190735 48.09265 0.095367 0.143051 For the 4th-order RK, the solution is t c k1 cmid k2 cmid k3 cend k4 0 10 2 20 1.5 17.5 1.625 26.25 1.1875 1.57291710 25.72917 1.213542 31.79688 0.910156 30.27995 0.986003 35.58919 0.72054 0.954420 35.27317 0.736342 38.95487 0.552256 38.03445 0.598278 41.25594 0.437203 0.57910230 41.06419 0.446791 43.29814 0.335093 42.73965 0.363017 44.69436 0.265282 0.35138240 44.57801 0.2711 45.93351 0.203325 45.59463 0.220268 46.78069 0.160965 0.21320850 46.71009 0.164495 47.53257 0.123371 47.32695 0.133652 48.04662 0.097669 0.129369 A plot of both solutions along with the analytical result is displayed below: 0255002040AnalyticalRK-4RK-2 28.2 The mass-balance equations can be written as 1130.16 0.06 1dcccdt 2120.2 0.2dcccdt 3230.05 0.25 4dcccdt 434 50.0875 0.125 0.0375dccc cdt 51250.04 0.02 0.06dccccdt Selected solution results (Euler’s method) are displayed below, along with a plot of the results. t c1 c2 c3 c4c5dc1/dt dc2/dt dc3/dt dc4/dt dc5/dt0 0.0000 0.0000 0.0000 0.0000 0.0000 1.0000 0.0000 4.0000 0.0000 0.00001 1.0000 0.0000 4.0000 0.0000 0.0000 1.0800 0.2000 3.0000 0.3500 0.04002 2.0800 0.2000 7.0000 0.3500 0.0400 1.0872 0.3760 2.2600 0.5703 0.08482 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. 3 3.1672 0.5760 9.2600 0.9203 0.1248 1.0488 0.5182 1.7138 0.6999 0.13074 4.2160 1.0942 10.9738 1.6201 0.2555 0.9839 0.6244 1.3113 0.7673 0.17525 5.1999 1.7186 12.2851 2.3874 0.4307 0.9051 0.6963 1.0147 0.7927 0.2165- - - 76 13.2416 13.2395 18.6473 16.8515 12.9430 0.0002 0.0004 0.0001 0.0106 0.017977 13.2418 13.2399 18.6475 16.8621 12.9609 0.0002 0.0004 0.0001 0.0099 0.016878 13.2419 13.2403 18.6476 16.8720 12.9777 0.0001 0.0003 0.0001 0.0093 0.015879 13.2421 13.2406 18.6477 16.8813 12.9935 0.0001 0.0003 0.0001 0.0088 0.014980 13.2422 13.2409 18.6478 16.8901 13.0084 0.0001 0.0003 0.0001 0.0082 0.0140051015200 1020304050607080c1 c2 c3 c4 c5 Finally, MATLAB can be used to determine the eigenvalues and eigenvectors: >> a=[.16 -.06 0 0 0;-.2 .2 0 0 0;0 -.05 .25 0 0;0 0 -.0875 .125 -.0375;-.04 -.02 0 0 .06] a = 0.1600 -0.0600 0 0 0 -0.2000 0.2000 0 0 0 0 -0.0500 0.2500 0 0 0 0 -0.0875 0.1250 -0.0375 -0.0400 -0.0200 0 0 0.0600 >> [v,d]=eig(a) v = 0 0 0 -0.1039 0.2604 0 0 0 -0.1582 -0.5701 0 0.8192 0 -0.0436 0.6892 1.0000 -0.5735 0.4997 0.4956 -0.3635 0 0 0.8662 0.8466 0.0043 d = 0.1250 0 0 0 0 0 0.2500 0 0 0 0 0 0.0600 0 0 0 0 0 0.0686 0 0 0 0 0 0.2914 28.3 Substituting the parameters into the differential equation gives 214.583333 0.0833333 0.15dcccdt The mid-point method can be applied with the result:3 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. The results are approaching a steady-state value of 9.586267. Challenge question: The steady state form (i.e., dc/dt = 0) of the equation is 0 = 14.58333 – 0.08333c – 0.15c2, which can be solved for 9.586267 and –10.1418. Thus, there is a negative root. If we put in the initial y value as 10.1418 (or higher precision), the solution will stay at the negative root for awhile until roundoff errors may lead to the solution going unstable. If we use an initial value less than the negative root, the solution will also go unstable. However, if we pick a value that is slightly higher (as per machine precision), it will gravitate towards the positive root. For example if we use 10.14 -15-10-50510150 5 10 15 20 This kind of behavior will also occur for higher initial conditions. For example, using an initial condition of 16 gives, 051015200 5 10 15 20 However, if we start to use even higher initial conditions, the solution will again go unstable. For example, if we use an initial condition of 17, the result is4 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. 010203040500 5 10 15 20 Therefore, it looks like initial conditions roughly in the range from 10.14 up to about 16.5 will yield stable solutions that converge on the steady-state solution of 9.586. Note that this range depends on our choice of step size. 28.4 The first steps of the solution are shown below along with a plot.
View Full Document
Unlocking...