Chapt_06 - 6 1. h=0.1 Numerical Solutions of Ordinary...

Info iconThis preview shows page 1. Sign up to view the full content.

View Full Document Right Arrow Icon
This is the end of the preview. Sign up to access the rest of the document.

Unformatted text preview: 6 1. h=0.1 Numerical Solutions of Ordinary Differential Equations EXERCISES 6.1 Euler Methods and Error Analysis h=0.05 2. yn 5.0000 4.4475 3.9763 3.5751 3.2342 2.9452 2.7009 2.4952 2.3226 2.1786 2.0592 h=0.1 h=0.05 xn 1.00 1.10 1.20 1.30 1.40 1.50 yn 5.0000 3.9900 3.2546 2.7236 2.3451 2.0801 xn 1.00 1.05 1.10 1.15 1.20 1.25 1.30 1.35 1.40 1.45 1.50 xn 0.00 0.10 0.20 0.30 0.40 0.50 yn 2.0000 1.6600 1.4172 1.2541 1.1564 1.1122 xn 0.00 0.05 0.10 0.15 0.20 0.25 0.30 0.35 0.40 0.45 0.50 yn 2.0000 1.8150 1.6571 1.5237 1.4124 1.3212 1.2482 1.1916 1.1499 1.1217 1.1056 3. h=0.1 h=0.05 4. yn 0.0000 0.0501 0.1004 0.1512 0.2028 0.2554 0.3095 0.3652 0.4230 0.4832 0.5465 h=0.1 h=0.05 xn 0.00 0.10 0.20 0.30 0.40 0.50 yn 0.0000 0.1005 0.2030 0.3098 0.4234 0.5470 xn 0.00 0.05 0.10 0.15 0.20 0.25 0.30 0.35 0.40 0.45 0.50 xn 0.00 0.10 0.20 0.30 0.40 0.50 yn 1.0000 1.1110 1.2515 1.4361 1.6880 2.0488 xn 0.00 0.05 0.10 0.15 0.20 0.25 0.30 0.35 0.40 0.45 0.50 yn 1.0000 1.0526 1.1113 1.1775 1.2526 1.3388 1.4387 1.5556 1.6939 1.8598 2.0619 317 6.1 5. Euler Methods and Error Analysis h=0.1 h=0.05 6. yn 0.0000 0.0488 0.0953 0.1397 0.1823 0.2231 0.2623 0.3001 0.3364 0.3715 0.4054 h=0.1 h=0.05 xn 0.00 0.10 0.20 0.30 0.40 0.50 yn 0.0000 0.0952 0.1822 0.2622 0.3363 0.4053 xn 0.00 0.05 0.10 0.15 0.20 0.25 0.30 0.35 0.40 0.45 0.50 xn 0.00 0.10 0.20 0.30 0.40 0.50 yn 0.0000 0.0050 0.0200 0.0451 0.0805 0.1266 xn 0.00 0.05 0.10 0.15 0.20 0.25 0.30 0.35 0.40 0.45 0.50 yn 0.0000 0.0013 0.0050 0.0113 0.0200 0.0313 0.0451 0.0615 0.0805 0.1022 0.1266 7. h=0.1 h=0.05 8. yn 0.5000 0.5116 0.5214 0.5294 0.5359 0.5408 0.5444 0.5469 0.5484 0.5492 0.5495 h=0.1 h=0.05 xn 0.00 0.10 0.20 0.30 0.40 0.50 yn 0.5000 0.5215 0.5362 0.5449 0.5490 0.5503 xn 0.00 0.05 0.10 0.15 0.20 0.25 0.30 0.35 0.40 0.45 0.50 xn 0.00 0.10 0.20 0.30 0.40 0.50 yn 1.0000 1.1079 1.2337 1.3806 1.5529 1.7557 xn 0.00 0.05 0.10 0.15 0.20 0.25 0.30 0.35 0.40 0.45 0.50 yn 1.0000 1.0519 1.1079 1.1684 1.2337 1.3043 1.3807 1.4634 1.5530 1.6503 1.7560 9. h=0.1 h=0.05 10. yn 1.0000 1.0024 1.0100 1.0228 1.0414 1.0663 1.0984 1.1389 1.1895 1.2526 1.3315 h=0.1 h=0.05 xn 1.00 1.10 1.20 1.30 1.40 1.50 yn 1.0000 1.0095 1.0404 1.0967 1.1866 1.3260 xn 1.00 1.05 1.10 1.15 1.20 1.25 1.30 1.35 1.40 1.45 1.50 xn 0.00 0.10 0.20 0.30 0.40 0.50 yn 0.5000 0.5250 0.5498 0.5744 0.5986 0.6224 xn 0.00 0.05 0.10 0.15 0.20 0.25 0.30 0.35 0.40 0.45 0.50 yn 0.5000 0.5125 0.5250 0.5374 0.5498 0.5622 0.5744 0.5866 0.5987 0.6106 0.6224 318 6.1 Euler Methods and Error Analysis 11. To obtain the analytic solution use the substitution u = x + y − 1. The resulting differential equation in u(x) will be separable. h=0.1 h=0.05 xn 0.00 0.10 0.20 0.30 0.40 0.50 yn 2.0000 2.1220 2.3049 2.5858 3.0378 3.8254 Actual Value 2.0000 2.1230 2.3085 2.5958 3.0650 3.9082 xn 0.00 0.05 0.10 0.15 0.20 0.25 0.30 0.35 0.40 0.45 0.50 yn 2.0000 2.0553 2.1228 2.2056 2.3075 2.4342 2.5931 2.7953 3.0574 3.4057 3.8840 Actual Value 2.0000 2.1230 2.3085 2.5958 3.0650 3.9082 2.5958 2.7997 3.0650 3.4189 3.9082 12. (a) y 20 15 10 5 1.1 1.2 1.3 1.4 x (b) xn 1.00 1.10 1.20 1.30 1.40 Euler 1.0000 1.2000 1.4938 1.9711 2.9060 Imp . Euler 1.0000 1.2469 1.6430 2.4042 4.5085 13. (a) Using Euler’s method we obtain y (0.1) ≈ y1 = 1.2. (b) Using y = 4e2x we see that the local truncation error is y (c) (0.1)2 h2 = 4e2c = 0.02e2c . 2 2 Since e2x is an increasing function, e2c ≤ e2(0.1) = e0.2 for 0 ≤ c ≤ 0.1. Thus an upper bound for the local truncation error is 0.02e0.2 = 0.0244. (c) Since y (0.1) = e0.2 = 1.2214, the actual error is y (0.1) − y1 = 0.0214, which is less than 0.0244. (d) Using Euler’s method with h = 0.05 we obtain y (0.1) ≈ y2 = 1.21. (e) The error in (d) is 1.2214 − 1.21 = 0.0114. With global truncation error O(h), when the step size is halved we expect the error for h = 0.05 to be one-half the error when h = 0.1. Comparing 0.0114 with 0.214 we see that this is the case. 14. (a) Using the improved Euler’s method we obtain y (0.1) ≈ y1 = 1.22. (b) Using y = 8e2x we see that the local truncation error is y (c) (0.1)3 h3 = 8e2c = 0.001333e2c . 6 6 Since e2x is an increasing function, e2c ≤ e2(0.1) = e0.2 for 0 ≤ c ≤ 0.1. Thus an upper bound for the local truncation error is 0.001333e0.2 = 0.001628. (c) Since y (0.1) = e0.2 = 1.221403, the actual error is y (0.1) − y1 = 0.001403 which is less than 0.001628. (d) Using the improved Euler’s method with h = 0.05 we obtain y (0.1) ≈ y2 = 1.221025. 319 6.1 Euler Methods and Error Analysis (e) The error in (d) is 1.221403 − 1.221025 = 0.000378. With global truncation error O(h2 ), when the step size is halved we expect the error for h = 0.05 to be one-fourth the error for h = 0.1. Comparing 0.000378 with 0.001403 we see that this is the case. 15. (a) Using Euler’s method we obtain y (0.1) ≈ y1 = 0.8. (b) Using y = 5e−2x we see that the local truncation error is 5e−2c (0.1)2 = 0.025e−2c . 2 Since e−2x is a decreasing function, e−2c ≤ e0 = 1 for 0 ≤ c ≤ 0.1. Thus an upper bound for the local truncation error is 0.025(1) = 0.025. (c) Since y (0.1) = 0.8234, the actual error is y (0.1) − y1 = 0.0234, which is less than 0.025. (d) Using Euler’s method with h = 0.05 we obtain y (0.1) ≈ y2 = 0.8125. (e) The error in (d) is 0.8234 − 0.8125 = 0.0109. With global truncation error O(h), when the step size is halved we expect the error for h = 0.05 to be one-half the error when h = 0.1. Comparing 0.0109 with 0.0234 we see that this is the case. 16. (a) Using the improved Euler’s method we obtain y (0.1) ≈ y1 = 0.825. (b) Using y = −10e−2x we see that the local truncation error is 10e−2c (0.1)3 = 0.001667e−2c . 6 Since e−2x is a decreasing function, e−2c ≤ e0 = 1 for 0 ≤ c ≤ 0.1. Thus an upper bound for the local truncation error is 0.001667(1) = 0.001667. (c) Since y (0.1) = 0.823413, the actual error is y (0.1) − y1 = 0.001587, which is less than 0.001667. (d) Using the improved Euler’s method with h = 0.05 we obtain y (0.1) ≈ y2 = 0.823781. (e) The error in (d) is |0.823413 − 0.8237181| = 0.000305. With global truncation error O(h2 ), when the step size is halved we expect the error for h = 0.05 to be one-fourth the error when h = 0.1. Comparing 0.000305 with 0.001587 we see that this is the case. 17. (a) Using y = 38e−3(x−1) we see that the local truncation error is y (c) h2 h2 = 38e−3(c−1) = 19h2 e−3(c−1) . 2 2 (b) Since e−3(x−1) is a decreasing function for 1 ≤ x ≤ 1.5, e−3(c−1) ≤ e−3(1−1) = 1 for 1 ≤ c ≤ 1.5 and y (c) h2 ≤ 19(0.1)2 (1) = 0.19. 2 (c) Using Euler’s method with h = 0.1 we obtain y (1.5) ≈ 1.8207. With h = 0.05 we obtain y (1.5) ≈ 1.9424. (d) Since y (1.5) = 2.0532, the error for h = 0.1 is E0.1 = 0.2325, while the error for h = 0.05 is E0.05 = 0.1109. With global truncation error O(h) we expect E0.1 /E0.05 ≈ 2. We actually have E0.1 /E0.05 = 2.10. 18. (a) Using y = −114e−3(x−1) we see that the local truncation error is y (c) h3 h3 = 114e−3(x−1) = 19h3 e−3(c−1) . 6 6 (b) Since e−3(x−1) is a decreasing function for 1 ≤ x ≤ 1.5, e−3(c−1) ≤ e−3(1−1) = 1 for 1 ≤ c ≤ 1.5 and y (c) h3 ≤ 19(0.1)3 (1) = 0.019. 6 320 6.2 Runge-Kutta Methods (c) Using the improved Euler’s method with h = 0.1 we obtain y (1.5) ≈ 2.080108. With h = 0.05 we obtain y (1.5) ≈ 2.059166. (d) Since y (1.5) = 2.053216, the error for h = 0.1 is E0.1 = 0.026892, while the error for h = 0.05 is E0.05 = 0.005950. With global truncation error O(h2 ) we expect E0.1 /E0.05 ≈ 4. We actually have E0.1 /E0.05 = 4.52. 19. (a) Using y = −1/(x + 1)2 we see that the local truncation error is y (c) h2 h2 1 = . 2 (c + 1)2 2 (b) Since 1/(x + 1)2 is a decreasing function for 0 ≤ x ≤ 0.5, 1/(c + 1)2 ≤ 1/(0 + 1)2 = 1 for 0 ≤ c ≤ 0.5 and y (c) (0.1)2 h2 ≤ (1) = 0.005. 2 2 (c) Using Euler’s method with h = 0.1 we obtain y (0.5) ≈ 0.4198. With h = 0.05 we obtain y (0.5) ≈ 0.4124. (d) Since y (0.5) = 0.4055, the error for h = 0.1 is E0.1 = 0.0143, while the error for h = 0.05 is E0.05 = 0.0069. With global truncation error O(h) we expect E0.1 /E0.05 ≈ 2. We actually have E0.1 /E0.05 = 2.06. 20. (a) Using y = 2/(x + 1)3 we see that the local truncation error is y (c) h3 1 h3 = . 6 (c + 1)3 3 (b) Since 1/(x + 1)3 is a decreasing function for 0 ≤ x ≤ 0.5, 1/(c + 1)3 ≤ 1/(0 + 1)3 = 1 for 0 ≤ c ≤ 0.5 and y (c) h3 (0.1)3 ≤ (1) = 0.000333. 6 3 (c) Using the improved Euler’s method with h = 0.1 we obtain y (0.5) ≈ 0.405281. With h = 0.05 we obtain y (0.5) ≈ 0.405419. (d) Since y (0.5) = 0.405465, the error for h = 0.1 is E0.1 = 0.000184, while the error for h = 0.05 is E0.05 = 0.000046. With global truncation error O(h2 ) we expect E0.1 /E0.05 ≈ 4. We actually have E0.1 /E0.05 = 3.98. ∗ ∗ 21. Because yn+1 depends on yn and is used to determine yn+1 , all of the yn cannot be computed at one time ∗ independently of the corresponding yn values. For example, the computation of y4 involves the value of y3 . EXERCISES 6.2 Runge-Kutta Methods 1. xn 0.00 0.10 0.20 0.30 0.40 0.50 yn 2.0000 2.1230 2.3085 2.5958 3.0649 3.9078 Actual Value 2.0000 2.1230 2.3085 2.5958 3.0650 3.9082 321 6.2 Runge-Kutta Methods 2. In this problem we use h = 0.1. Substituting w2 = in (4) in the text, we obtain w1 = 1 − w2 = 1 , 4 α= 2 1 =, 2w2 3 3 4 into the equations xn and β = 2 1 =. 2w2 3 The resulting second-order Runge-Kutta method is 1 3 k1 + k2 4 4 0.00 0.10 0.20 0.30 0.40 0.50 Second Order Improved Runge Kutta Euler 2.0000 2.0000 2.1213 2.1220 2.3030 2.3049 2.5814 2.5858 3.0277 3.0378 3.8002 3.8254 yn+1 = yn + h where = yn + h (k1 + 3k2 ) 4 k1 = f (xn , yn ) k2 = f 2 2 xn + h, yn + hk1 . 3 3 The table compares the values obtained using this second-order Runge-Kutta method with the values obtained using the improved Euler’s method. 3. xn 1.00 1.10 1.20 1.30 1.40 1.50 yn 5.0000 3.9724 3.2284 2.6945 2.3163 2.0533 4. xn 0.00 0.10 0.20 0.30 0.40 0.50 yn 2.0000 1.6562 1.4110 1.2465 1.1480 1.1037 5. xn 0.00 0.10 0.20 0.30 0.40 0.50 yn 0.0000 0.1003 0.2027 0.3093 0.4228 0.5463 6. xn 0.00 0.10 0.20 0.30 0.40 0.50 yn 1.0000 1.1115 1.2530 1.4397 1.6961 2.0670 7. xn 0.00 0.10 0.20 0.30 0.40 0.50 yn 0.0000 0.0953 0.1823 0.2624 0.3365 0.4055 8. xn 0.00 0.10 0.20 0.30 0.40 0.50 yn 0.0000 0.0050 0.0200 0.0451 0.0805 0.1266 9. xn 0.00 0.10 0.20 0.30 0.40 0.50 yn 0.5000 0.5213 0.5358 0.5443 0.5482 0.5493 10. xn 0.00 0.10 0.20 0.30 0.40 0.50 yn 1.0000 1.1079 1.2337 1.3807 1.5531 1.7561 322 6.2 11. xn 1.00 1.10 1.20 1.30 1.40 1.50 yn 1.0000 1.0101 1.0417 1.0989 1.1905 1.3333 Runge-Kutta Methods 12. xn 0.00 0.10 0.20 0.30 0.40 0.50 yn 0.5000 0.5250 0.5498 0.5744 0.5987 0.6225 13. (a) Write the equation in the form dv = 32 − 0.125v 2 = f (t, v ). dt (b) v 40 30 20 10 t tn 0.0 1.0 2.0 3.0 4.0 5.0 vn 0.0000 25.2570 32.9390 34.9770 35.5500 35.7130 1 2 3 4 5 6 (c) Separating variables and using partial fractions we have 1 √ 2 32 and √ 1 √ √ 32 − 1 √ 0.125 v +√ 32 + 1 √ 0.125 v dv = dt 2 32 0.125 Since v (0) = 0 we find c = 0. Solving for v we obtain √ √ 16 5 (e 3.2 t − 1) √ v (t) = e 3 .2 t + 1 √ √ √ √ ln | 32 + 0.125 v | − ln | 32 − 0.125 v | = t + c. and v (5) ≈ 35.7678. Alternatively, the solution can be expressed as v (t) = mg tanh k kg t. m 14. (a) t days A observed A approximated A(t) 1 2.78 1.93 2 13.53 12.50 3 36.30 36.46 4 47.50 47.23 5 49.40 49.00 (b) From the graph we estimate A(1) ≈ 1.68, A(2) ≈ 13.2, A(3) ≈ 36.8, A(4) ≈ 46.9, and A(5) ≈ 48.9. 50 40 30 20 10 0 1 2 3 4 5 t 323 6.2 Runge-Kutta Methods (c) Let α = 2.128 and β = 0.0432. Separating variables we obtain dA = dt A(α − βA) 1 α 1 β + A α − βA dA = dt 1 [ln A − ln(α − βA)] = t + c α A ln = α(t + c) α − βA A = eα(t+c) α − βA A = αeα(t+c) − βAeα(t+c) 1 + βeα(t+c) A = αeα(t+c) . Thus A(t) = αeα(t+c) α α = = . β + e−αc e−αt 1 + βeα(t+c) β + e−α(t+c) From A(0) = 0.24 we obtain 0.24 = so that e−αc = α/0.24 − β ≈ 8.8235 and A(t) ≈ t days A observed A actual 1 2.78 1.93 2 13.53 12.50 α β + e−αc 2.128 . 0.0432 + 8.8235e−2.128t 4 47.50 47.23 5 49.40 49.00 3 36.30 36.46 15. (a) xn 1.00 1.05 1.10 1.15 1.20 1.25 1.30 1.35 1.40 h 0.05 1.0000 1.1112 1.2511 1.4348 1.6934 2.1047 2.9560 7.8981 1.0608 1015 h 0.1 1.0000 1.2511 1.6934 2.9425 (b) y 20 15 10 5 1.1 1.2 1.3 1.4 x 903.0282 16. (a) Using the RK4 method we obtain y (0.1) ≈ y1 = 1.2214. (b) Using y (5) (x) = 32e2x we see that the local truncation error is y (5) (c) (0.1)5 h5 = 32e2c = 0.000002667e2c . 120 120 324 6.2 Runge-Kutta Methods Since e2x is an increasing function, e2c ≤ e2(0.1) = e0.2 for 0 ≤ c ≤ 0.1. Thus an upper bound for the local truncation error is 0.000002667e0.2 = 0.000003257. (c) Since y (0.1) = e0.2 = 1.221402758, the actual error is y (0.1) − y1 = 0.000002758 which is less than 0.000003257. (d) Using the RK4 formula with h = 0.05 we obtain y (0.1) ≈ y2 = 1.221402571. (e) The error in (d) is 1.221402758 − 1.221402571 = 0.000000187. With global truncation error O(h4 ), when the step size is halved we expect the error for h = 0.05 to be one-sixteenth the error for h = 0.1. Comparing 0.000000187 with 0.000002758 we see that this is the case. 17. (a) Using the RK4 method we obtain y (0.1) ≈ y1 = 0.823416667. (b) Using y (5) (x) = −40e−2x we see that the local truncation error is 40e−2c (0.1)5 = 0.000003333. 120 Since e−2x is a decreasing function, e−2c ≤ e0 = 1 for 0 ≤ c ≤ 0.1. Thus an upper bound for the local truncation error is 0.000003333(1) = 0.000003333. (c) Since y (0.1) = 0.823413441, the actual error is |y (0.1) − y1 | = 0.000003225, which is less than 0.000003333. (d) Using the RK4 method with h = 0.05 we obtain y (0.1) ≈ y2 = 0.823413627. (e) The error in (d) is |0.823413441 − 0.823413627| = 0.000000185. With global truncation error O(h4 ), when the step size is halved we expect the error for h = 0.05 to be one-sixteenth the error when h = 0.1. Comparing 0.000000185 with 0.000003225 we see that this is the case. 18. (a) Using y (5) = −1026e−3(x−1) we see that the local truncation error is y (5) (c) h5 = 8.55h5 e−3(c−1) . 120 (b) Since e−3(x−1) is a decreasing function for 1 ≤ x ≤ 1.5, e−3(c−1) ≤ e−3(1−1) = 1 for 1 ≤ c ≤ 1.5 and y (5) (c) h5 ≤ 8.55(0.1)5 (1) = 0.0000855. 120 With h = 0.05 we obtain (c) Using the RK4 method with h = 0.1 we obtain y (1.5) ≈ 2.053338827. y (1.5) ≈ 2.053222989. 19. (a) Using y (5) = 24/(x + 1)5 we see that the local truncation error is y (5) (c) h5 h5 1 = . 55 120 (c + 1) (b) Since 1/(x + 1)5 is a decreasing function for 0 ≤ x ≤ 0.5, 1/(c + 1)5 ≤ 1/(0 + 1)5 = 1 for 0 ≤ c ≤ 0.5 and y (5) (c) (0.1)5 h5 ≤ (1) = 0.000002. 5 5 With h = 0.05 we obtain (c) Using the RK4 method with h = 0.1 we obtain y (0.5) ≈ 0.405465168. y (0.5) ≈ 0.405465111. 20. Each step of Euler’s method requires only 1 function evaluation, while each step of the improved Euler’s method ∗ requires 2 function evaluations – once at (xn , yn ) and again at (xn+1 , yn+1 ). The second-order Runge-Kutta methods require 2 function evaluations per step, while the RK4 method requires 4 function evaluations per step. To compare the methods we approximate the solution of y = (x + y − 1)2 , y (0) = 2, at x = 0.2 using h = 0.1 325 6.2 Runge-Kutta Methods for the Runge-Kutta method, h = 0.05 for the improved Euler’s method, and h = 0.025 for Euler’s method. For each method a total of 8 function evaluations is required. By comparing with the exact solution we see that the RK4 method appears to still give the most accurate result. xn 0.000 0.025 0.050 0.075 0.100 0.125 0.150 0.175 0.200 Imp . Euler Euler h 0.025 h 0.05 2.0000 2.0000 2.0250 2.0526 2.0553 2.0830 2.1165 2.1228 2.1535 2.1943 2.2056 2.2395 2.2895 2.3075 RK4 h 0.1 2.0000 Actual 2.0000 2.0263 2.0554 2.0875 2.1230 2.1624 2.2061 2.2546 2.3085 y 2.1230 2.3085 21. (a) For y + y = 10 sin 3x an integrating factor is ex so that dx [e y ] = 10ex sin 3x =⇒ ex y = ex sin 3x − 3ex cos 3x + c dx =⇒ y = sin 3x − 3 cos 3x + ce−x . When x = 0, y = 0, so 0 = −3 + c and c = 3. The solution is y = sin 3x − 3 cos 3x + 3e−x . Using Newton’s method we find that x = 1.53235 is the only positive root in [0, 2]. 5 2 x −5 (b) Using the RK4 method with h = 0.1 we obtain the table of values shown. These values are used to obtain an interpolating function in Mathematica. The graph of the interpolating function is shown. Using Mathematica’s root finding capability we see that the only positive root in [0, 2] is x = 1.53236. xn yn xn yn 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 0.0000 0.1440 0.5448 1.1409 1.8559 2.6049 3.3019 3.8675 4.2356 4.3593 4.2147 1.0 1.1 1.2 1.3 1.4 1.5 1.6 1.7 1.8 1.9 2.0 4.2147 3.8033 3.1513 2.3076 1.3390 0.3243 0.6530 1.5117 2.1809 2.6061 2.7539 y 5 2 x −5 326 6.3 Multistep Methods EXERCISES 6.3 Multistep Methods In the tables in this section “ABM” stands for Adams-Bashforth-Moulton. 1. Writing the differential equation in the form y − y = x − 1 we see that an integrating factor is e− so that d −x [e y ] = (x − 1)e−x dx and y = ex (−xe−x + c) = −x + cex . dx = e−x , From y (0) = 1 we find c = 1, so the solution of the initial-value problem is y = −x + ex . Actual values of the analytic solution above are compared with the approximated values in the table. xn 0.0 0.2 0.4 0.6 0.8 yn 1.00000000 1.02140000 1.09181796 1.22210646 1.42552788 Actual 1.00000000 1.02140276 1.09182470 1.22211880 1.42554093 init. cond. RK4 RK4 RK4 ABM 2. The following program is written in Mathematica. It uses the Adams-Bashforth-Moulton method to approximate the solution of the initial-value problem y = x + y − 1, y (0) = 1, on the interval [0, 1]. Clear[f, x, y, h, a, b, y0]; f[x , y ]:= x + y - 1; h = 0.2; a = 0; y0 = 1; b = 1; f[x, y] Clear[k1, k2, k3, k4, x, y, u, v] x = u[0] = a; y = v[0] = y0; n = 0; While[x < a + 3h, n = n + 1; (* define the differential equation *) (* set the step size *) (* set the initial condition and the interval *) (* display the DE *) (* use RK4 to compute the first 3 values after y(0) *) k1 = f[x, y]; k2 = f[x + h/2, y + h k1/2]; k3 = f[x + h/2, y + h k2/2]; k4 = f[x + h, y + h k3]; x = x + h; y = y + (h/6)(k1 + 2k2 + 2k3 + k4); u[n] = x; v[n] = y]; 327 6.3 Multistep Methods While[x ≤ b, (* use Adams-Bashforth-Moulton *) p3 = f[u[n - 3], v[n - 3]]; p2 = f[u[n - 2], v[n - 2]]; p1 = f[u[n - 1], v[n - 1]]; p0 = f[u[n], v[n]]; pred = y + (h/24)(55p0 - 59p1 + 37p2 - 9p3); (* predictor *) x = x + h; p4 = f[x, pred]; y = y + (h/24)(9p4 + 19p0 - 5p1 + p2); n = n + 1; u[n] = x; v[n] = y] (*display the table *) TableForm[Prepend[Table[{u[n], v[n]}, {n, 0, (b-a)/h}], {"x(n)", "y(n)"}]]; (* corrector *) ∗ 3. The first predictor is y4 = 0.73318477. xn 0.0 0.2 0.4 0.6 0.8 yn 1.00000000 0.73280000 0.64608032 0.65851653 0.72319464 init. cond. RK4 RK4 RK4 ABM ∗ 4. The first predictor is y4 = 1.21092217. xn 0.0 0.2 0.4 0.6 0.8 yn 2.00000000 1.41120000 1.14830848 1.10390600 1.20486982 init. cond. RK4 RK4 RK4 ABM ∗ 5. The first predictor for h = 0.2 is y4 = 1.02343488. xn 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 h 0.2 0.00000000 init. cond. 0.20270741 RK4 0.42278899 RK4 0.68413340 RK4 1.02969040 ABM 1.55685960 ABM h 0.1 0.00000000 0.10033459 0.20270988 0.30933604 0.42279808 0.54631491 0.68416105 0.84233188 1.02971420 1.26028800 1.55762558 init. cond. RK4 RK4 RK4 ABM ABM ABM ABM ABM ABM ABM 328 6.3 ∗ 6. The first predictor for h = 0.2 is y4 = 3.34828434. Multistep Methods xn 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 h 0.2 1.00000000 init. cond. 1.44139950 RK4 1.97190167 RK4 2.60280694 RK4 3.34860927 ABM 4.22797875 ABM h 0.1 1.00000000 1.21017082 1.44140511 1.69487942 1.97191536 2.27400341 2.60283209 2.96031780 3.34863769 3.77026548 4.22801028 init. cond. RK4 RK4 RK4 ABM ABM ABM ABM ABM ABM ABM ∗ 7. The first predictor for h = 0.2 is y4 = 0.13618654. xn 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 h 0.2 0.00000000 init. cond. 0.00262739 RK4 0.02005764 RK4 0.06296284 RK4 0.13598600 ABM 0.23854783 ABM h 0.1 0.00000000 0.00033209 0.00262486 0.00868768 0.02004821 0.03787884 0.06294717 0.09563116 0.13596515 0.18370712 0.23841344 init. cond. RK4 RK4 RK4 ABM ABM ABM ABM ABM ABM ABM ∗ 8. The first predictor for h = 0.2 is y4 = 2.61796154. xn 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 h 0.2 1.00000000 init. cond. 1.23369623 RK4 1.55308554 RK4 1.99610329 RK4 2.62136177 ABM 3.52079042 ABM h 0.1 1.00000000 1.10793839 1.23369772 1.38068454 1.55309381 1.75610064 1.99612995 2.28119129 2.62131818 3.02914333 3.52065536 init. cond. RK4 RK4 RK4 ABM ABM ABM ABM ABM ABM ABM 329 6.3 Multistep Methods 6.4 Higher-Order Equations and Systems EXERCISES 6.4 Higher-Order Equations and Systems 1. The substitution y = u leads to the iteration formulas yn+1 = yn + hun , The initial conditions are y0 = −2 and u0 = 1. Then y1 = y0 + 0.1u0 = −2 + 0.1(1) = −1.9 u1 = u0 + 0.1(4u0 − 4y0 ) = 1 + 0.1(4 + 8) = 2.2 y2 = y1 + 0.1u1 = −1.9 + 0.1(2.2) = −1.68. The general solution of the differential equation is y = c1 e2x + c2 xe2x . From the initial conditions we find c1 = −2 and c2 = 5. Thus y = −2e2x + 5xe2x and y (0.2) ≈ 1.4918. 2. The substitution y = u leads to the iteration formulas yn+1 = yn + hun , un+1 = un + h 2 2 u n − 2 yn . x x un+1 = un + h(4un − 4yn ). The initial conditions are y0 = 4 and u0 = 9. Then y1 = y0 + 0.1u0 = 4 + 0.1(9) = 4.9 u1 = u0 + 0.1 2 2 u 0 − y0 1 1 = 9 + 0.1[2(9) − 2(4)] = 10 y2 = y1 + 0.1u1 = 4.9 + 0.1(10) = 5.9. The general solution of the Cauchy-Euler differential equation is y = c1 x + c2 x2 . From the initial conditions we find c1 = −1 and c2 = 5. Thus y = −x + 5x2 and y (1.2) = 6. 3. The substitution y = u leads to the system y = u, u = 4u − 4y. xn 0.0 0.1 0.2 h 0.2 yn 2.0000 1.4928 h 0.2 un 1.0000 4.4731 h 0.1 yn 2.0000 1.8321 1.4919 h 0.1 un 1.0000 2.4427 4.4753 Using formula (4) in the text with x corresponding to t, y corresponding to x, and u corresponding to y , we obtain the table shown. 4. The substitution y = u leads to the system y = u, 2 2 u = u − 2 y. x x xn 1.0 1.1 1.2 h 0.2 yn 4.0000 6.0001 h 0.2 un 9.0000 11.0002 h 0.1 yn 4.0000 4.9500 6.0000 h 0.1 un 9.0000 10.0000 11.0000 Using formula (4) in the text with x corresponding to t, y corresponding to x, and u corresponding to y , we obtain the table shown. 330 6.4 Higher-Order Equations and Systems 5. The substitution y = u leads to the system y = u, u = 2u − 2y + et cos t. xn 0.0 0.1 0.2 h 0.2 yn 1.0000 1.4640 h 0.2 un 2.0000 2.6594 h 0.1 yn 1.0000 1.2155 1.4640 h 0.1 un 2.0000 2.3150 2.6594 Using formula (4) in the text with y corresponding to x and u corresponding to y , we obtain the table shown. 6. Using h = 0.1, the RK4 method for a system, and a numerical solver, we obtain i1 7 6 5 4 3 2 1 1 2 3 4 5t i2 7 6 5 4 3 2 1 1 2 3 4 5t tn 0.0 0.1 0.2 0.3 0.4 0.5 h 0.2 i 1n 0.0000 2.5000 2.8125 2.0703 0.6104 1.5619 h 0.2 i 3n 0.0000 3.7500 5.7813 7.4023 9.1919 11.4877 7. tn 0.0 0.1 0.2 h 0.2 xn 6.0000 8.3055 h 0.2 yn 2.0000 3.4199 h 0.1 xn 6.0000 7.0731 8.3055 h 0.1 yn 2.0000 2.6524 3.4199 x,y 20 xHtL 15 10 yHtL 0.5 1 1.5 2t 5 8. x,y tn 0.0 0.1 0.2 h 0.2 xn 1.0000 2.0785 h 0.2 yn 1.0000 3.3382 h 0.1 xn 1.0000 1.4006 2.0845 h 0.1 yn 1.0000 1.8963 3.3502 50 yHtL 40 30 20 xHtL 10 0.5 1 1.5 2t 331 6.4 9. Higher-Order Equations and Systems x,y tn 0.0 0.1 0.2 h 0.2 xn 3.0000 3.9123 h 0.2 yn 5.0000 4.2857 h 0.1 xn 3.0000 3.4790 3.9123 h 0.1 yn 5.0000 4.6707 4.2857 30 25 20 15 10 5 xHtL 5 -5 10 15 yHtL 20 25 30 t 10. tn 0.0 0.1 0.2 h 0.2 xn 0.5000 2.1589 h 0.2 yn 0.2000 2.3279 h 0.1 xn 0.5000 1.0207 2.1904 h 0.1 yn 0.2000 1.0115 2.3592 x,y 2 xHtL 1.5 1 yHtL 0.1 0.2 0.5 t 0.3 11. Solving for x and y we obtain the system x = −2x + y + 5t y = 2x + y − 2t. x,y 1 2 3 xHtL 4t -5 tn 0.0 0.1 0.2 h 0.2 xn 1.0000 0.4179 h 0.2 yn 2.0000 2.1824 h 0.1 xn 1.0000 0.6594 0.4173 h 0.1 yn 2.0000 2.0476 2.1821 -10 -15 -20 yHtL 332 6.5 Second-Order Boundary-Value Problems 12. Solving for x and y we obtain the system x,y 1 x = y − 3t2 + 2t − 5 2 1 y = − y + 3t2 + 2t + 5. 2 h 0.2 xn 3.0000 1.9867 h 0.2 yn 1.0000 0.0933 h 0.1 xn 3.0000 2.4727 1.9867 h 0.1 yn 1.0000 0.4527 0.0933 60 40 20 yHtL 2 -20 -40 -60 4 6 8t xHtL tn 0.0 0.1 0.2 EXERCISES 6.5 Second-Order Boundary-Value Problems 1. We identify P (x) = 0, Q(x) = 9, f (x) = 0, and h = (2 − 0)/4 = 0.5. Then the finite difference equation is yi+1 + 0.25yi + yi−1 = 0. The solution of the corresponding linear system gives x y 0.0 0.5 1.0 1.5 4.0000 -5.6774 -2.5807 6.3226 2.0 1.0000 2. We identify P (x) = 0, Q(x) = −1, f (x) = x2 , and h = (1 − 0)/4 = 0.25. Then the finite difference equation is yi+1 − 2.0625yi + yi−1 = 0.0625x2 . i The solution of the corresponding linear system gives x y 0.00 0.25 0.50 0.75 1.00 0.0000 -0.0172 -0.0316 -0.0324 0.0000 3. We identify P (x) = 2, Q(x) = 1, f (x) = 5x, and h = (1 − 0)/5 = 0.2. Then the finite difference equation is 1.2yi+1 − 1.96yi + 0.8yi−1 = 0.04(5xi ). The solution of the corresponding linear system gives x y 0.0 0.2 0.4 0.6 0.8 1.0 0.0000 -0.2259 -0.3356 -0.3308 -0.2167 0.0000 333 6.5 Second-Order Boundary-Value Problems 4. We identify P (x) = −10, Q(x) = 25, f (x) = 1, and h = (1 − 0)/5 = 0.2. Then the finite difference equation is −yi + 2yi−1 = 0.04. The solution of the corresponding linear system gives x y 0.0 1.0000 0.2 1.9600 0.4 3.8800 0.6 0.8 1.0 7.7200 15.4000 0.0000 5. We identify P (x) = −4, Q(x) = 4, f (x) = (1 + x)e2x , and h = (1 − 0)/6 = 0.1667. Then the finite difference equation is 0.6667yi+1 − 1.8889yi + 1.3333yi−1 = 0.2778(1 + xi )e2xi . The solution of the corresponding linear system gives x y 0.0000 3.0000 0.1667 3.3751 0.3333 3.6306 0.5000 3.6448 0.6667 3.2355 0.8333 2.1411 1.0000 0.0000 √ 6. We identify P (x) = 5, Q(x) = 0, f (x) = 4 x , and h = (2 − 1)/6 = 0.1667. Then the finite difference equation is √ 1.4167yi+1 − 2yi + 0.5833yi−1 = 0.2778(4 xi ). The solution of the corresponding linear system gives x y 1.0000 1.1667 1.3333 1.5000 1.6667 1.8333 2.0000 1.0000 -0.5918 -1.1626 -1.3070 -1.2704 -1.1541 -1.0000 7. We identify P (x) = 3/x, Q(x) = 3/x2 , f (x) = 0, and h = (2 − 1)/8 = 0.125. Then the finite difference equation is 0.1875 0.0469 0.1875 1+ yi+1 + −2 + yi + 1 − yi−1 = 0. xi x2 xi i The solution of the corresponding linear system gives x y 1.000 5.0000 1.125 3.8842 1.250 2.9640 1.375 2.2064 1.500 1.5826 1.625 1.0681 1.750 0.6430 1.875 0.2913 2.000 0.0000 8. We identify P (x) = −1/x, Q(x) = x−2 , f (x) = ln x/x2 , and h = (2 − 1)/8 = 0.125. Then the finite difference equation is 0.0625 0.0156 0.0625 1− yi+1 + −2 + yi + 1 + yi−1 = 0.0156 ln xi . 2 xi xi xi The solution of the corresponding linear system gives x y 1.000 1.125 1.250 1.375 1.500 1.625 1.750 1.875 2.000 0.0000 -0.1988 -0.4168 -0.6510 -0.8992 -1.1594 -1.4304 -1.7109 -2.0000 9. We identify P (x) = 1 − x, Q(x) = x, f (x) = x, and h = (1 − 0)/10 = 0.1. Then the finite difference equation is [1 + 0.05(1 − xi )]yi+1 + [−2 + 0.01xi ]yi + [1 − 0.05(1 − xi )]yi−1 = 0.01xi . The solution of the corresponding linear system gives x y 0.0 0.0000 0.1 0.2660 0.2 0.5097 0.3 0.7357 0.4 0.9471 0.5 0.6 1.1465 1.3353 0.7 0.8 1.5149 1.6855 0.9 1.8474 1.0 2.0000 10. We identify P (x) = x, Q(x) = 1, f (x) = x, and h = (1 − 0)/10 = 0.1. Then the finite difference equation is (1 + 0.05xi )yi+1 − 1.99yi + (1 − 0.05xi )yi−1 = 0.01xi . 334 6.5 The solution of the corresponding linear system gives x y 0.0 1.0000 0.1 0.8929 0.2 0.7789 0.3 0.6615 0.4 0.5440 0.5 0.6 0.4296 0.3216 0.7 0.8 0.2225 0.1347 Second-Order Boundary-Value Problems 0.9 0.0601 1.0 0.0000 11. We identify P (x) = 0, Q(x) = −4, f (x) = 0, and h = (1 − 0)/8 = 0.125. Then the finite difference equation is yi+1 − 2.0625yi + yi−1 = 0. The solution of the corresponding linear system gives x y 0.000 0.0000 0.125 0.3492 0.250 0.7202 0.375 1.1363 0.500 1.6233 0.625 2.2118 0.750 2.9386 0.875 3.8490 1.000 5.0000 12. We identify P (r) = 2/r, Q(r) = 0, f (r) = 0, and h = (4 − 1)/6 = 0.5. Then the finite difference equation is 1+ 0.5 ri ui+1 − 2ui + 1 − 0.5 ri ui−1 = 0. The solution of the corresponding linear system gives r 1.0 1.5 2.0 2.5 3.0 3.5 4.0 u 50.0000 72.2222 83.3333 90.0000 94.4444 97.6190 100.0000 13. (a) The difference equation 1+ h h Pi yi+1 + (−2 + h2 Qi )yi + 1 − Pi yi−1 = h2 fi 2 2 is the same as equation (8) in the text. The equations are the same because the derivation was based only on the differential equation, not the boundary conditions. If we allow i to range from 0 to n − 1 we obtain n equations in the n + 1 unknowns y−1 , y0 , y1 , . . . , yn−1 . Since yn is one of the given boundary conditions, it is not an unknown. (b) Identifying y0 = y (0), y−1 = y (0 − h), and y1 = y (0 + h) we have from equation (5) in the text 1 [y1 − y−1 ] = y (0) = 1 2h The difference equation corresponding to i = 0, 1+ h h P0 y1 + (−2 + h2 Q0 )y0 + 1 − P0 y−1 = h2 f0 2 2 or y1 − y−1 = 2h. becomes, with y−1 = y1 − 2h, 1+ or 2y1 + (−2 + h2 Q0 )y0 = h2 f0 + 2h − P0 . Alternatively, we may simply add the equation y1 − y−1 = 2h to the list of n difference equations obtaining n + 1 equations in the n + 1 unknowns y−1 , y0 , y1 , . . . , yn−1 . (c) Using n = 5 we obtain x 0.0 0.2 0.4 0.6 0.8 1.0 y -2.2755 -2.0755 -1.8589 -1.6126 -1.3275 -1.0000 h h P0 y1 + (−2 + h2 Q0 )y0 + 1 − P0 (y1 − 2h) = h2 f0 2 2 335 6.5 Second-Order Boundary-Value Problems 14. Using h = 0.1 and, after shooting a few times, y (0) = 0.43535 we obtain the following table with the RK4 method. x y 0.0 1.00000 0.1 1.04561 0.2 1.09492 0.3 1.14714 0.4 1.20131 0.7 1.36392 0.5 1.25633 0.8 1.41388 0.6 1.31096 0.9 1.45962 1.0 1.50003 CHAPTER 6 REVIEW EXERCISES 1. xn 1.00 1.05 1.10 1.15 1.20 1.25 1.30 1.35 1.40 1.45 1.50 Euler h 0.1 2.0000 2.1386 2.3097 2.5136 2.7504 3.0201 Euler Imp . Euler Imp . Euler h 0.05 h 0.1 h 0.05 2.0000 2.0000 2.0000 2.0693 2.0735 2.1469 2.1549 2.1554 2.2328 2.2459 2.3272 2.3439 2.3450 2.4299 2.4527 2.5409 2.5672 2.5689 2.6604 2.6937 2.7883 2.8246 2.8269 2.9245 2.9686 3.0690 3.1157 3.1187 RK4 h 0.1 2.0000 2.1556 2.3454 2.5695 2.8278 3.1197 RK4 h 0.05 2.0000 2.0736 2.1556 2.2462 2.3454 2.4532 2.5695 2.6944 2.8278 2.9696 3.1197 2. xn 0.00 0.05 0.10 0.15 0.20 0.25 0.30 0.35 0.40 0.45 0.50 Euler h 0.1 0.0000 0.1000 0.2010 0.3049 0.4135 0.5279 Euler Imp . Euler Imp . Euler h 0.05 h 0.1 h 0.05 0.0000 0.0000 0.0000 0.0500 0.0501 0.1001 0.1005 0.1004 0.1506 0.1512 0.2017 0.2030 0.2027 0.2537 0.2552 0.3067 0.3092 0.3088 0.3610 0.3638 0.4167 0.4207 0.4202 0.4739 0.4782 0.5327 0.5382 0.5378 RK4 h 0.1 0.0000 0.1003 0.2026 0.3087 0.4201 0.5376 RK4 h 0.05 0.0000 0.0500 0.1003 0.1511 0.2026 0.2551 0.3087 0.3637 0.4201 0.4781 0.5376 336 CHAPTER 6 REVIEW EXERCISES 3. xn 0.50 0.55 0.60 0.65 0.70 0.75 0.80 0.85 0.90 0.95 1.00 Euler h 0.1 0.5000 0.6000 0.7095 0.8283 0.9559 1.0921 Imp . Euler Imp . Euler Euler h 0.05 h 0.1 h 0.05 0.5000 0.5000 0.5000 0.5500 0.5512 0.6024 0.6048 0.6049 0.6573 0.6609 0.7144 0.7191 0.7193 0.7739 0.7800 0.8356 0.8427 0.8430 0.8996 0.9082 0.9657 0.9752 0.9755 1.0340 1.0451 1.1044 1.1163 1.1168 RK4 h 0.1 0.5000 0.6049 0.7194 0.8431 0.9757 1.1169 RK4 h 0.05 0.5000 0.5512 0.6049 0.6610 0.7194 0.7801 0.8431 0.9083 0.9757 1.0452 1.1169 4. xn 1.00 1.05 1.10 1.15 1.20 1.25 1.30 1.35 1.40 1.45 1.50 Euler h 0.1 1.0000 1.2000 1.4760 1.8710 2.4643 3.4165 Euler Imp . Euler Imp . Euler h 0.05 h 0.1 h 0.05 1.0000 1.0000 1.0000 1.1000 1.1091 1.2183 1.2380 1.2405 1.3595 1.4010 1.5300 1.5910 1.6001 1.7389 1.8523 1.9988 2.1524 2.1799 2.3284 2.6197 2.7567 3.1458 3.2360 3.3296 4.1528 4.1253 5.2510 5.6404 RK4 h 0.1 1.0000 1.2415 1.6036 2.1909 3.2745 5.8338 RK4 h 0.05 1.0000 1.1095 1.2415 1.4029 1.6036 1.8586 2.1911 2.6401 3.2755 4.2363 5.8446 5. Using yn+1 = yn + hun , un+1 = un + h(2xn + 1)yn , y1 = y0 + 0.1u0 = 3 + (0.1)1 = 3.1 u1 = u0 + 0.1(2x0 + 1)y0 = 1 + 0.1(1)3 = 1.3 y2 = y1 + 0.1u1 = 3.1 + 0.1(1.3) = 3.23. ∗ 6. The first predictor is y3 = 1.14822731. y0 = 3 u0 = 1 we obtain (when h = 0.2) y1 = y (0.2) = y0 + hu0 = 3 + (0.2)1 = 3.2. When h = 0.1 we have xn 0.0 0.1 0.2 0.3 0.4 yn 2.00000000 1.65620000 1.41097281 1.24645047 1.14796764 init. cond. RK4 RK4 RK4 ABM 7. Using x0 = 1, y0 = 2, and h = 0.1 we have x1 = x0 + h(x0 + y0 ) = 1 + 0.1(1 + 2) = 1.3 y1 = y0 + h(x0 − y0 ) = 2 + 0.1(1 − 2) = 1.9 and 337 CHAPTER 6 REVIEW EXERCISES x2 = x1 + h(x1 + y1 ) = 1.3 + 0.1(1.3 + 1.9) = 1.62 y2 = y1 + h(x1 − y1 ) = 1.9 + 0.1(1.3 − 1.9) = 1.84. Thus, x(0.2) ≈ 1.62 and y (0.2) ≈ 1.84. 8. We identify P (x) = 0, Q(x) = 6.55(1 + x), f (x) = 1, and h = (1 − 0)/10 = 0.1. Then the finite difference equation is yi+1 + [−2 + 0.0655(1 + xi )]yi + yi−1 = 0.001 or yi+1 + (0.0655xi − 1.9345)yi + yi−1 = 0.001. The solution of the corresponding linear system gives x y 0.0 0.0000 0.1 4.1987 0.2 0.3 0.4 0.5 0.6 8.1049 11.3840 13.7038 14.7770 14.4083 0.7 12.5396 0.8 9.2847 0.9 4.9450 1.0 0.0000 338 ...
View Full Document

This note was uploaded on 02/17/2010 for the course MATHEMATIC MAS201 taught by Professor Xingqin during the Spring '10 term at Korea Advanced Institute of Science and Technology.

Ask a homework question - tutors are online