Bài giảng Phương pháp tính - Phương trình vi phân - Đậu Thế Phiệt

PHƯƠNG TRÌNH VI PHÂN Bài giảng điện tử Ngày 6 tháng 12 năm 2016 PHƯƠNG TRÌNH VI PHÂN Ngày 6 tháng 12 năm 2016 1 / 54 CuuDuongThanCong.com https://fb.com/tailieudientucntt Ta xét bài toán cơ bản về dao động của con lắc đơn C H A P T E R 5 Initial-Value Problems for Ordinary Differential Equations Introduction The motion of a swinging pendulum under certain simplifying assumptions is described by the second-order differential equation d2θ dt2 + g L sin θ = 0, L θ where L is the length of the pendulum, g ≈ 32.17 ft/s2 is the gravitational constant of the earth, and θ is the angle the pendulum makes with the vertical. If, in addition, we specify the position of the pendulum when the motion begins, θ(t0) = θ0, and its velocity at that point, θ ′(t0) = θ ′0, we have what is called an initial-value problem. For small values of θ , the approximation θ ≈ sin θ can be used to simplify this problem to the linear initial-value problem d2θ dt2 + g L θ = 0, θ(t0) = θ0, θ ′(t0) = θ ′0. This problem can be solved by a standard differential-equation technique. For larger values of θ , the assumption that θ = sin θ is not reasonable so approximation methods must be used. A problem of this type is considered in Exercise 8 of Section 5.9. Any textbook on ordinary differential equations details a number of methods for ex- plicitly finding solutions to first-order initial-value problems. In practice, however, few of the problems originating from the study of physical phenomena can be solved exactly. 259 Copyright 2010 Cengage Learning. All Rights Reserved. May not be copied, scanned, or duplicated, in whole or in part. Due to electronic rights, some third party content may be suppressed from the eBook and/or eChapter(s). Editorial review has deemed that any suppressed content does not materially affect the overall learning experience. Cengage Learning reserves the right to remove additional content at any time if subsequent rights restrictions require it. xác định bởi phương trình vi phân bậc hai d2θ dt2 + g L sin θ = 0 với L là chiều dài con lắc, g là hằng số hấp dẫn của trái đất, θ là góc tạo bởi con lắc và trục thẳng đứng. Ta xét vị trí ban đầu của con lắc khi bắt đầu dao động là θ(t0) = θ0 và vận tốc ban đầu tại điểm này là θ′(t0) = θ′0, ta có bài toán giá trị đầu. PHƯƠNG TRÌNH VI PHÂN Ngày 6 tháng 12 năm 2016 2 / 54 CuuDuongThanCong.com https://fb.com/tailieudientucntt Với giá trị θ nhỏ, ta xấp xỉ θ ≈ sin θ, khi đó bài toán trở thành tuyến tính d2θ dt2 + g L θ = 0, θ(t0) = θ0, θ ′(t0) = θ′0 Bài toán này có thể giải bằng các phương pháp quen thuộc. Tuy nhiên với giá trị θ lớn, ta không thể giả thiết θ = sin θ. Để tìm nghiệm cho bài toán này, ta cần sử dụng các phương pháp xấp xỉ nghiệm. PHƯƠNG TRÌNH VI PHÂN Ngày 6 tháng 12 năm 2016 3 / 54 CuuDuongThanCong.com https://fb.com/tailieudientucntt Bài toán Cauchy Đặt vấn đề Bài toán Cauchy Ta xét bài toán giá trị đầu bậc nhất, bài toán Cauchy,{ y ′(t) = f (t, y(t)), a 6 t 6 b, y(a) = α (1) với y = y(t) là hàm cần tìm, khả vi trên đoạn [a, b], y0 là giá trị ban đầu cho trước của y(t) tại t = a. PHƯƠNG TRÌNH VI PHÂN Ngày 6 tháng 12 năm 2016 4 / 54 CuuDuongThanCong.com https://fb.com/tailieudientucntt Bài toán Cauchy Đặt vấn đề Đối với bài toán Cauchy (1) ta chỉ có thể tìm được nghiệm đúng của một số phương trình đơn giản, còn đối với trường hợp f (x , y) có dạng bất kỳ thì nói chung không có phương pháp giải. Ngoài ra, trong những trường hợp có thể tìm ra nghiệm đúng của bài toán Cauchy (1) quá phức tạp thì người ta cũng ít dùng. Vì vậy, việc tìm những phương pháp giải gần đúng bài toán Cauchy có vai trò rất quan trọng trong thực tế. PHƯƠNG TRÌNH VI PHÂN Ngày 6 tháng 12 năm 2016 5 / 54 CuuDuongThanCong.com https://fb.com/tailieudientucntt Bài toán Cauchy Công thức Euler Công thức Euler Để tìm nghiệm gần đúng của bài toán (1) ta chia đoạn [a, b] thành n đoạn nhỏ bằng nhau với h = b − a n . Khi đó các điểm nút là t0 = a, tk = t0 + kh, k = 0, 1, 2, . . . , n, tn = b. Giả sử y(t) là nghiệm duy nhất của bài toán (1), có đạo hàm đến cấp 2 liên tục trên đoạn [a, b]. Khi đó với mỗi k = 0, 1, 2, . . . , n − 1 theo công thức khai triển Taylor trên đoạn [tk , tk+1], ta có y(tk+1) = y(tk) + y ′(tk)(tk+1 − tk) + y ′′(ξk)(tk+1 − tk) 2 2 , với ξk ∈ (tk , tk+1). PHƯƠNG TRÌNH VI PHÂN Ngày 6 tháng 12 năm 2016 6 / 54 CuuDuongThanCong.com https://fb.com/tailieudientucntt Bài toán Cauchy Công thức Euler Vì y = y(t) là nghiệm của phương trình (1) và h = tk+1 − tk nên ta có y(tk+1) = y(tk) + h.f (tk , yk) + h2 2 y ′′(ξk) Bằng cách bỏ đi phần dư, ta xấp xỉ yk ≈ y(tk) với k = 1, 2, . . . n, ta có công thức Euler y0 = α yk+1 ≈ yk + hf (tk , yk), với k = 0, 1, 2, . . . , n − 1. PHƯƠNG TRÌNH VI PHÂN Ngày 6 tháng 12 năm 2016 7 / 54 CuuDuongThanCong.com https://fb.com/tailieudientucntt Bài toán Cauchy Công thức Euler Ý nghĩa hình học của phương pháp Euler 268 C H A P T E R 5 Initial-Value Problems for Ordinary Differential Equations The graph of the function highlighting y(ti) is shown in Figure 5.2. One step in Euler’s method appears in Figure 5.3, and a series of steps appears in Figure 5.4. Figure 5.2 t y y(tN) ! y(b) y" ! f (t, y), y(a) ! α y(t2) y(t1) y(t0) ! α t0 ! a t1 t2 tN ! b. . . . . . Figure 5.3 w1 Slope y"(a) ! f (a, α) y t y" ! f (t, y), y(a) ! α α t0 ! a t1 t2 tN ! b. . . Figure 5.4 w1 y t α t 0 ! a t1 t2 tN ! b y(b) w2 wN y" ! f (t, y), y(a) ! α . . . Example 1 Euler’s method was used in the first illustration with h = 0.5 to approximate the solution to the initial-value problem y′ = y − t2 + 1, 0 ≤ t ≤ 2, y(0) = 0.5. Use Algorithm 5.1 with N = 10 to determine approximations, and compare these with the exact values given by y(t) = (t + 1)2 − 0.5et . Solution With N = 10 we have h = 0.2, ti = 0.2i, w0 = 0.5, and wi+1 = wi + h(wi − t2i + 1) = wi + 0.2[wi − 0.04i2 + 1] = 1.2wi − 0.008i2 + 0.2, for i = 0, 1, . . . , 9. So w1 = 1.2(0.5)− 0.008(0)2 + 0.2 = 0.8; w2 = 1.2(0.8)− 0.008(1)2 + 0.2 = 1.152; and so on. Table 5.1 shows the comparison between the approximate values at ti and the actual values. Copyright 2010 Cengage Learning. All Rights Reserved. May not be copied, scanned, or duplicated, in whole or in part. Due to electronic rights, some third party content may be suppressed from the eBook and/or eChapter(s). Editorial review has deemed that any suppressed content does not materially affect the overall learning experience. Cengage Learning reserves the right to remove additional content at any time if subsequent rights restrictions require it. 268 C H A P T E R 5 Initial-Value Problems for Ordinary Differential Equations The graph of the function highlighting y(ti) is shown in Figure 5.2. One step in Euler’s method appears in Figure 5.3, and a series of steps appears in Figure 5.4. Figure 5.2 t y y(tN) ! y(b) y" ! f (t, y), y(a) ! α y(t2) y(t1) y(t0) ! α t0 ! a t1 t2 tN ! b. . . . . . Figure 5.3 w1 Slope y"(a) ! f (a, α) y t y" ! f (t, y), y(a) ! α α t0 ! a t1 t2 tN ! b. . . Figure 5.4 w1 y t α t 0 ! a t1 t2 tN ! b y(b) w2 wN y" ! f (t, y), y(a) ! α . . . Example 1 Euler’s method was used in the first illustration with h = 0.5 to approximate the solution to the initial-value problem y′ = y − t2 + 1, 0 ≤ t ≤ 2, y(0) = 0.5. Use Algorithm 5.1 with N = 10 to determine approximations, and compare these with the exact values given by y(t) = (t + 1)2 − 0.5et . Solution With N = 10 we have h = 0.2, ti = 0.2i, w0 = 0.5, and wi+1 = wi + h(wi − t2i + 1) = wi + 0.2[wi − 0.04i2 1] = 1.2wi − 0.008i2 + 0.2, for i = 0, 1, . . . , 9. So w1 = 1.2(0.5)− 0.008(0)2 + 0.2 = 0.8; w2 = 1.2(0.8)− 0.008(1)2 + 0.2 = 1.152; and so on. Table 5.1 shows the comparison between the approximate values at ti and the actual values. Copyright 2010 Cengage Learning. All Rights Reserved. May not be copied, scanned, or duplicated, in whole or in part. Due to electronic rights, some third party content may be suppressed from the eBook and/or eChapter(s). Editorial review has deemed that any suppressed content does not materially affect the overall learning experience. Cengage Learning reserves the right to remove additional content at any time if subsequent rights restrictions require it. Từ (t0, y0) = (a, α) thuộc đường cong y = y(t), kẻ tiếp tuyến với đường cong (có hệ số góc là y ′(a) = f (a, α)). Đường tiếp tuyến sẽ cắt t = t1 tại y1 chính là giá trị gần đúng của y(t1). Tại (t1, y1), ta kẻ đường thẳng với hệ số góc f (t1, y1) cắt t = t2 tại y2 là giá trị gần đúng của y(t2). PHƯƠNG TRÌNH VI PHÂN Ngày 6 tháng 12 năm 2016 8 / 54 CuuDuongThanCong.com https://fb.com/tailieudientucntt Bài toán Cauchy Công thức Euler Ví dụ Sử dụng phương pháp Euler để xấp xỉ nghiệm của bài toán Cauchy{ y ′(x) = y − t2 + 1, 0 6 t 6 2, y(0) = 0.5 với n = 10. Tại những điểm nút chia so sánh giá trị gần đúng với giá trị chính xác, biết nghiệm chính xác của bài toán là y(t) = (t + 1)2 − 0.5et . PHƯƠNG TRÌNH VI PHÂN Ngày 6 tháng 12 năm 2016 9 / 54 CuuDuongThanCong.com https://fb.com/tailieudientucntt Bài toán Cauchy Công thức Euler Giải. Với n = 10 thì h = 2− 0 10 = 0.2, tk = 0.2k, y0 = 0.5. Công thức tính nghiệm gần đúng là yk+1 = yk + h(yk − t2k + 1) với k = 0, 1, . . . , 9. Bấm máy. Y = Y + 0.2(Y − X 2 + 1) : X = X + 0.2 1 CALC Y = 0.5 =,X = 0 = 2 Y =,X = 0.2 = PHƯƠNG TRÌNH VI PHÂN Ngày 6 tháng 12 năm 2016 10 / 54 CuuDuongThanCong.com https://fb.com/tailieudientucntt Bài toán Cauchy Công thức Euler Giải. Với n = 10 thì h = 2− 0 10 = 0.2, tk = 0.2k, y0 = 0.5. Công thức tính nghiệm gần đúng là yk+1 = yk + h(yk − t2k + 1) với k = 0, 1, . . . , 9. Bấm máy. Y = Y + 0.2(Y − X 2 + 1) : X = X + 0.2 1 CALC Y = 0.5 =,X = 0 = 2 Y =,X = 0.2 = PHƯƠNG TRÌNH VI PHÂN Ngày 6 tháng 12 năm 2016 10 / 54 CuuDuongThanCong.com https://fb.com/tailieudientucntt Bài toán Cauchy Công thức Euler k tk yk y(tk) |y(tk)− yk | 0 0.0 0.5000000 0.5000000 0.0000000 1 0.2 0.8000000 0.8292986 0.0292986 2 0.4 1.1520000 1.2140877 0.0620877 3 0.6 1.5504000 1.6489406 0.0985406 4 0.8 1.9884800 2.1272295 0.1387495 5 1.0 2.4581760 2.6408591 0.1826831 6 1.2 2.9498112 3.1799415 0.2301303 7 1.4 3.4517734 3.7324000 0.2806266 8 1.6 3.9501281 4.2834838 0.3333557 9 1.8 4.4281538 4.8151763 0.3870225 10 2.0 4.8657845 5.3054720 0.4396874 PHƯƠNG TRÌNH VI PHÂN Ngày 6 tháng 12 năm 2016 11 / 54 CuuDuongThanCong.com https://fb.com/tailieudientucntt Bài toán Cauchy Công thức Euler Sai số của công thức Euler Giả sử f là hàm liên tục và thỏa điều kiện |f (t, y1)− f (t, y2)| ≤ L|y1 − y2| với hằng số L > 0, và tồn tại M thỏa y ′′(t) ≤ M với t ∈ [a, b]. Khi đó với y(t) là nghiệm chính xác của bài toán giá trị đầu y ′(t) = f (t, y), a ≤ t ≤ b, y(a) = α và y0, y1, . . . , yn là nghiệm xấp xỉ của bài toán cho bởi công thức Euler, khi đó với mỗi k = 0, 1, . . . , n |y(tk)− yk | ≤ hM 2L [eL(tk−a) − 1] PHƯƠNG TRÌNH VI PHÂN Ngày 6 tháng 12 năm 2016 12 / 54 CuuDuongThanCong.com https://fb.com/tailieudientucntt Bài toán Cauchy Công thức Euler cải tiến Công thức Euler cải tiến Trong công thức Euler, thay f (tk , yk) bởi f (tk , yk) + f (tk+1, yk+1) 2 ta được công thức Euler cải tiến y(tk+1) ≈ yk+1 = yk + h f (tk , yk) + f (tk+1, yk+1) 2 , với k = 0, 1, 2, . . . , n − 1. Việc tính toán theo công thức Euler cải tiến rất phức tạp vì cả 2 vế đều chứa yk+1 là ẩn cần tìm. Để đơn giản ta thay yk+1 ở vế phải bởi yk + hf (tk , yk). PHƯƠNG TRÌNH VI PHÂN Ngày 6 tháng 12 năm 2016 13 / 54 CuuDuongThanCong.com https://fb.com/tailieudientucntt Bài toán Cauchy Công thức Euler cải tiến Lúc này ta có công thức y(xk+1) ≈ yk+1 = yk + h f (tk , yk) + f (tk+1, yk + hf (tk , yk)) 2 , k = 0, 1, 2, . . . , n − 1. PHƯƠNG TRÌNH VI PHÂN Ngày 6 tháng 12 năm 2016 14 / 54 CuuDuongThanCong.com https://fb.com/tailieudientucntt Bài toán Cauchy Công thức Euler cải tiến Ví dụ Sử dụng phương pháp Euler cải tiến để xấp xỉ nghiệm của bài toán Cauchy{ y ′(t) = y − t2 + 1, 0 6 t 6 2, y(0) = 0.5 với n = 10. Tại những điểm nút chia so sánh giá trị gần đúng với giá trị chính xác, biết nghiệm chính xác của bài toán là y(t) = (t + 1)2 − 0.5et . PHƯƠNG TRÌNH VI PHÂN Ngày 6 tháng 12 năm 2016 15 / 54 CuuDuongThanCong.com https://fb.com/tailieudientucntt Bài toán Cauchy Công thức Euler cải tiến Với n = 10 thì h = 2− 0 10 = 0.2, y0 = 0.5. Công thức tính nghiệm gần đúng là yk+1 = yk + h f (tk , yk) + f (tk+1, yk + hf (tk , yk)) 2 với k = 0, 1, . . . , 9. Bấm máy. Y = Y + 0.1× (Y − X 2 + 1 + Y + 0.2(Y − X 2 + 1)− (X + 0.2)2 + 1) : X = X + 0.2 1 CALC Y = 0.5 = X = 0 = 2 Y =,X = 0.2 PHƯƠNG TRÌNH VI PHÂN Ngày 6 tháng 12 năm 2016 16 / 54 CuuDuongThanCong.com https://fb.com/tailieudientucntt Bài toán Cauchy Công thức Euler cải tiến k tk yk y(tk) |y(tk)− yk | 0 0.0 0.5 0.5000000 0.0000000 1 0.2 0.826 0.8292986 0.0032986 2 0.4 1.20692 1.2140877 0.0071677 3 0.6 1.6372424 1.6489406 0.0116982 4 0.8 2.110235728 2.1272295 0.0169938 5 1.0 2.617687588 2.6408591 0.0231715 6 1.2 3.149578858 3.1799415 0.0303627 7 1.4 3.693686206 3.7324000 0.0387138 8 1.6 4.235097172 4.2834838 0.0483866 9 1.8 4.755618549 4.8151763 0.0595577 10 2.0 5.23305463 5.3054720 0.0724173 PHƯƠNG TRÌNH VI PHÂN Ngày 6 tháng 12 năm 2016 17 / 54 CuuDuongThanCong.com https://fb.com/tailieudientucntt Bài toán Cauchy Công thức Euler cải tiến Bài tập Sử dụng công thức Euler và công thức Euler cải tiến để xấp xỉ nghiệm của các bài toán sau 1 y ′ = te3t − 2y , 0 ≤ t ≤ 1, y(0) = 0 với h = 0.5. 2 y ′ = cos 2t + sin 3t, 0 ≤ t ≤ 1, y(0) = 1 với h = 0.25 3 y ′ = 1 + y/t, 1 ≤ t ≤ 2, y(1) = 2 với h = 0.25. 4 y ′ = 1 + t 1 + y , 1 ≤ t ≤ 2, y(1) = 2, với h = 0.25 PHƯƠNG TRÌNH VI PHÂN Ngày 6 tháng 12 năm 2016 18 / 54 CuuDuongThanCong.com https://fb.com/tailieudientucntt Bài toán Cauchy Công thức Runge-Kutta Công thức Runge- Kutta bậc hai Xét khai triển Taylor bậc hai của y(t), ta có y(tk+1) = y(tk) + hy ′(tk) + h2 2 y ′′(tk) + h3 3! y ′′′(ξ) = y(tk) + hf (tk , y(tk)) + h2 2 f ′(tk , y(tk)) + h3 3! y ′′′(ξ) Ta lại có f ′(tk , y(tk) = ∂f ∂t (tk , y(tk) + ∂f ∂y (tk , y(tk)).y ′(tk) và y ′(tk) = f (tk , y(tk)) y(tk+1) ≈ y(tk)+h { f (tk , y(tk)) + h 2 ∂f ∂t (tk , y(tk)) + h 2 ∂f ∂y (tk , y(tk)).f (tk , y(tk)) } PHƯƠNG TRÌNH VI PHÂN Ngày 6 tháng 12 năm 2016 19 / 54 CuuDuongThanCong.com https://fb.com/tailieudientucntt Bài toán Cauchy Công thức Runge-Kutta Công thức Runge- Kutta bậc hai Xét khai triển Taylor bậc hai của y(t), ta có y(tk+1) = y(tk) + hy ′(tk) + h2 2 y ′′(tk) + h3 3! y ′′′(ξ) = y(tk) + hf (tk , y(tk)) + h2 2 f ′(tk , y(tk)) + h3 3! y ′′′(ξ) Ta lại có f ′(tk , y(tk) = ∂f ∂t (tk , y(tk) + ∂f ∂y (tk , y(tk)).y ′(tk) và y ′(tk) = f (tk , y(tk)) y(tk+1) ≈ y(tk)+h { f (tk , y(tk)) + h 2 ∂f ∂t (tk , y(tk)) + h 2 ∂f ∂y (tk , y(tk)).f (tk , y(tk)) } PHƯƠNG TRÌNH VI PHÂN Ngày 6 tháng 12 năm 2016 19 / 54 CuuDuongThanCong.com https://fb.com/tailieudientucntt Bài toán Cauchy Công thức Runge-Kutta Từ công thức khai triển Taylor a1f (tk + α, y(tk) + β) ≈ a1 [ f (tk , y(tk) + α ∂f ∂t (tk , y(tk) + β ∂f ∂y (tk , y(tk)) ] = a1f (tk , y(tk)) + a1α ∂f ∂t (tk , y(tk)) + a1β ∂f ∂y (tk , y(tk)) Chọn a1, α, β thỏa 1 = a1, h 2 = a1α, h 2 f (tk , y(tk)) = a1β ta được a1 = 1, α = h 2 , β = h 2 f (tk , y(tk)) Vậy y(tk+1) ≈ y(tk) + hf ( tk + h 2 , y(tk) + h 2 f (tk , y(tk) ) PHƯƠNG TRÌNH VI PHÂN Ngày 6 tháng 12 năm 2016 20 / 54 CuuDuongThanCong.com https://fb.com/tailieudientucntt Bài toán Cauchy Công thức Runge-Kutta Công thức điểm giữa y0 = α yk+1 = yk + h [ f ( tk + h 2 , yk + h 2 f (tk , yk) )] với k = 0, 1, . . . , n − 1 và sai số O(h2). PHƯƠNG TRÌNH VI PHÂN Ngày 6 tháng 12 năm 2016 21 / 54 CuuDuongThanCong.com https://fb.com/tailieudientucntt Bài toán Cauchy Công thức Runge-Kutta Nếu ta sử dụng a1f (t, y) + a2f (t + α, y + βf (t, y)) và đồng nhất các hệ số, ta có a1 = a2 = 1 2 và α = β = h. Ta có công thức Euler cải tiến Công thức Euler cải tiến y0 = α yk+1 = yk + h 2 [f (tk , yk) + f (tk+1, yk + hf (tk , yk))] với k = 0, 1, . . . , n − 1 và sai số O(h2). PHƯƠNG TRÌNH VI PHÂN Ngày 6 tháng 12 năm 2016 22 / 54 CuuDuongThanCong.com https://fb.com/tailieudientucntt Bài toán Cauchy Công thức Runge-Kutta Công thức Runge-Kutta có độ chính xác cao hơn công thức Euler, vì dùng khai triển Taylor nghiệm y = y(x) của bài toán (1) với nhiều số hạng hơn. Sử dụng quá trình xây dựng trên đối với công thức Taylor bâc cao hơn, ta có thể xây dựng Phương pháp Runge - Kutta với các bậc cao, và phổ biến nhất là bậc 4  yk+1 = y(tk + h) ≈ yk + 1 6 (K k1 + 2K k 2 + 2K k 3 + K k 4 ) K k1 = hf (tk , yk) K k2 = hf ( tk + h 2 , yk + K k1 2 ) K k3 = hf ( tk + h 2 , yk + K k2 2 ) K k4 = hf (tk + h, yk + K k 3 ) PHƯƠNG TRÌNH VI PHÂN Ngày 6 tháng 12 năm 2016 23 / 54 CuuDuongThanCong.com https://fb.com/tailieudientucntt Bài toán Cauchy Công thức Runge-Kutta Ví dụ Sử dụng phương pháp Runge-Kutta bậc 4 để xấp xỉ nghiệm của bài toán Cauchy { y ′(x) = y − t2 + 1, 0 6 t 6 2, y(0) = 0.5 với n = 10. Tại những điểm nút chia so sánh giá trị gần đúng với giá trị chính xác, biết nghiệm chính xác của bài toán là y(t) = (t + 1)2 − 0.5et . PHƯƠNG TRÌNH VI PHÂN Ngày 6 tháng 12 năm 2016 24 / 54 CuuDuongThanCong.com https://fb.com/tailieudientucntt Bài toán Cauchy Công thức Runge-Kutta Giải. Với n = 10 thì h = 2− 0 10 = 0.2, tk = 0.2k, y0 = 0.5. Ta có K k1 = hf (tk , yk) = 0.2(yk − 0.04k2 + 1), K k2 = hf ( tk + h 2 , yk + K k1 2 ) = h [ yk + K k1 2 − ( tk + h 2 )2 + 1 ] , K k3 = hf ( tk + h 2 , yk + K k2 2 ) = h [ yk + K k2 2 − ( tk + h 2 )2 + 1 ] , K k4 = hf (tk + h, yk + K k 3 ) = h [ yk + K k 3 − (xk + h)2 + 1 ] . Công thức tính nghiệm gần đúng là yk+1 = yk + 1 6 (K k1 + 2K k 2 + 2K k 3 + K k 4 ) với k = 0, 1, . . . , 9. PHƯƠNG TRÌNH VI PHÂN Ngày 6 tháng 12 năm 2016 25 / 54 CuuDuongThanCong.com https://fb.com/tailieudientucntt Bài toán Cauchy Công thức Runge-Kutta Bấm máy. 0.2(Y − X 2 + 1) Tính K 01 . CALC X = 0,Y = 0.5 ⇒ K 01 Shift-STO-A Tính K 02 . CALC X = 0 + 0.2÷ 2,Y = 0.5 + A÷ 2. ⇒ K 02 Shift-STO-B Tính K 03 . CALC X = 0 + 0.2÷ 2,Y = 0.5 + B ÷ 2. ⇒ K 03 Shift-STO-C Tính K 04 . CALC X = 0 + 0.2,Y = 0.5 + C . ⇒ K 04 Shift-STO-D PHƯƠNG TRÌNH VI PHÂN Ngày 6 tháng 12 năm 2016 26 / 54 CuuDuongThanCong.com https://fb.com/tailieudientucntt Bài toán Cauchy Công thức Runge-Kutta y(0.2) ≈ y1 = y0 + 1 6 (K 01 + 2K 0 2 + 2K 0 3 + K 0 4 ) = = 0.5 + 1 6 (A+ 2B + 2C + D) ≈ 0.8292933 Shift - STO - F PHƯƠNG TRÌNH VI PHÂN Ngày 6 tháng 12 năm 2016 27 / 54 CuuDuongThanCong.com https://fb.com/tailieudientucntt Bài toán Cauchy Công thức Runge-Kutta k tk yk y(tk) |y(tk)− yk | 0 0.0 0.5000000 0.5000000 0.0000000 1 0.2 0.8292933 0.8292986 0.0000053 2 0.4 1.2140762 1.2140877 0.0000114 3 0.6 1.6489220 1.6489406 0.0000186 4 0.8 2.1272027 2.1272295 0.0000269 5 1.0 2.6408227 2.6408591 0.0000364 6 1.2 3.1798942 3.1799415 0.0000474 7 1.4 3.7323401 3.7324000 0.0000599 8 1.6 4.2834095 4.2834838 0.0000743 9 1.8 44.8150857 4.8151763 0.0000906 10 2.0 5.3053630 5.3054720 0.0001089 PHƯƠNG TRÌNH VI PHÂN Ngày 6 tháng 12 năm 2016 28 / 54 CuuDuongThanCong.com https://fb.com/tailieudientucntt Bài toán Cauchy Công thức Runge-Kutta Bài tập Sử dụng công thức Runge - Kutta bậc 4 để xấp xỉ nghiệm cho các phương trình sau 1 y ′ = y t − (y t )2 với 1 ≤ t ≤ 2, y(1) = 1 và h = 0.1, nghiệm chính xác y(t) = t/(t + ln t). 2 y ′ = 1+ y t + (y t )2 với 1 ≤ t ≤ 3, y(1) = 0 và h = 0.2, nghiệm chính xác y(t) = t tan(ln t). 3 y ′ = −5y + 5t2 + 2t với 0 ≤ t ≤ 1, y(1) = 1/3 và h = 0.1, nghiệm chính xác y(t) = t2 + 1 3 e−5t . PHƯƠNG TRÌNH VI PHÂN Ngày 6 tháng 12 năm 2016 29 / 54 CuuDuongThanCong.com https://fb.com/tailieudientucntt Hệ phương trình vi phân Hệ phương trình vi phân Xét hệ phương trình vi phân x ′(t) = f (t, x(t), y(t)) y ′(t) = g(t, x(t), y(t)) x(t0) = α y(t0) = β t ∈ [t0, t0 + H] Chia đoạn [t0, t0 + H] thành n đoạn nhỏ bằng nhau có độ dài h = H n . Các điểm chia là tk = t0 + kh, k = 0, 1, . . . , n. Giá trị gần đúng tại điểm tk của x(t) là xk = x(tk), của y(t) là yk = y(tk) PHƯƠNG TRÌNH VI PHÂN Ngày 6 tháng 12 năm 2016 30 / 54 CuuDuongThanCong.com https://fb.com/tailieudientucntt Hệ phương trình vi phân Công thức Euler Công thức Euler  x(tk) ≈ xk = xk−1 + hf (tk−1, xk−1, yk−1) y(tk) ≈ yk = yk−1 + hg(tk−1, xk−1, yk−1) k = 1, 2, . . . , n PHƯƠNG TRÌNH VI PHÂN Ngày 6 tháng 12 năm 2016 31 / 54 CuuDuongThanCong.com https://fb.com/tailieudientucntt Hệ phương trình vi phân Công thức Euler cải tiến Công thức Euler cải tiến  K1x = hf (tk−1, xk−1, yk−1) K1y = hg(tk−1, xk−1, yk−1) K2x = hf (tk−1 + h, xk−1 + K1x , yk−1 + K1y ) K2y = hg(tk−1 + h, xk−1 + K1x , yk−1 + K1y ) x(tk) ≈ xk = xk−1 + 1 2 (K1x + K2x) y(tk) ≈ yk = yk−1 + 1 2 (K1y + K2y ) k = 1, 2, . . . , n PHƯƠNG TRÌNH VI PHÂN Ngày 6 tháng 12 năm 2016 32 / 54 CuuDuongThanCong.com https://fb.com/tailieudientucntt Hệ phương trình vi phân Công thức Runge-Kutta bậc bốn Công thức Runge-Kutta bậc bốn  K1x = hf (tk−1, xk−1, yk−1) K1y = hg(tk−1, xk−1, yk−1) K2x = hf ( tk−1 + h 2 , xk−1 + K1x 2 , yk−1 + K1y 2 ) K2y = hg ( tk−1 + h 2 , xk−1 + K1x 2 , yk−1 + K1y 2 ) K3x = hf ( tk−1 + h 2 , xk−1 + K2x 2 , yk−1 + K2y 2 ) K3y = hg ( tk−1 + h 2 , xk−1 + K2x 2 , yk−1 + K2y 2 ) K4x = hf (tk−1 + h, xk−1 + K3x , yk−1 + K3y ), K4y = hg(tk−1 + h, xk−1 + K3x , yk−1 + K3y ) x(tk) ≈ xk = xk−1 + 1 6 (K1x + 2K2x + 2K3x + K4x) y(tk) ≈ yk = yk−1 + 1 6 (K1y + 2K2y + 2K3y + K4y ) với k = 1, 2, . . . , n PHƯƠNG TRÌNH VI PHÂN Ngày 6 tháng 12 năm 2016 33 / 54 CuuDuongThanCong.com https://fb.com/tailieudientucntt Hệ phương trình vi phân Công thức Runge-Kutta bậc bốn Ví dụ Cho hệ  x ′(t) = tx − 2y + 1 y ′(t) = 2x + ty + sin t x(1) = 0.25 y(1) = 0.75 , t > 1 Sử dụng công thức Euler cải tiến để xấp xỉ giá trị của x(t) và y(t) tại t = 1.2 với bước h = 0.2. PHƯƠNG TRÌNH VI PHÂN Ngày 6 tháng 12 năm 2016 34 / 54 CuuDuongThanCong.com https://fb.com/tailieudientucntt Hệ phương trình vi phân Công thức Runge-Kutta bậc bốn Ta có t0 = 1, x0 = 0.25, y0 = 0.75, h = 0.2, f (t, x , y) = tx − 2y + 1, g(t, x , y) = 2x + ty + sin t K1x = hf (t0, x0, y0) = h(t0x0 − 2y0 + 1) = −0.05 K1y = hg(t0, x0, y0) = h(2x0 + t0y0 + sin t0) = 0.4183 K2x = hf (t0 + h, x0 + K1x , y0 + K1y ) = h[(t0 + h)(x0 + K1x)− 2(y0 + K1y ) + 1] = −0.2193 K2y = hg(t0 + h, x0 + K1x , y0 + K1y ) = h[2(x0 + K1x) + (t0 + h)(y0 + K1y ) + sin(t0 + h)] = 0.5468 x(1.2) ≈ x1 = x0 + 1 2 (K1x + K2x) = 0.1154 y(1.2) ≈ y1 = y0 + 1 2 (K1y + K2y ) = 1.2326 PHƯƠNG TRÌNH VI PHÂN Ngày 6 tháng 12 năm 2016 35 / 54 CuuDuongThanCong.com https://fb.com/tailieudientucntt Hệ phương trình vi phân Công thức Runge-Kutta bậc bốn Đối với phương trình vi phân cấp 2, x ′′(t) = F (t, x(t), x ′(t)) x(t0) = α x ′(t0) = β t ∈ [t0, t0 + H] được chuyển về hệ phương trình vi phân cấp 1 bằng cách đặt y(t) = x ′(t) x ′(t) = f (t, x(t), y(t)) = y y ′(t) = g(t, x(t), y(t)) = F (t, x , y) x(t0) = α y(t0) = β PHƯƠNG TRÌNH VI PHÂN Ngày 6 tháng 12 năm 2016 36 / 54 CuuDuongThanCong.com https://fb.com/tailieudientucntt Hệ phương trình vi phân Công thức Runge-Kutta bậc bốn Ví dụ Cho phương trình vi phân cấp 2 x ′′ − 2x ′ + 2x = e2t . sin t, (0 < t < 1) với điều kiện ban đầu x(0) = −0.4, x ′(0) = −0.6. Dùng phương pháp Runge-Kutta bậc 4 để xấp xỉ nghiệm gần đúng của phương trình với bước h = 0.1. So sánh kết quả thu được với nghiệm chính xác x(t) = 0.2e2t(sin t − 2 cos t), y(t) = x ′(t) = 0.2e2t(4 sin t − 3 cos t). PHƯƠNG TRÌNH VI PHÂN Ngày 6 tháng 12 năm 2016 37 / 54 CuuDuongThanCong.com https://fb.com/tailieudientucntt Hệ phương trình vi phân Công thức Runge-Kutta bậc bốn Đặt y(t) = x ′(t). Phương trình đã cho được biến đổi thành hệ x ′(t) = f (t, x(t), y(t)) = y y ′(t) = g(t, x(t), y(t)) = −2x + 2y + e2t sin t x(0) = −0.4 y(0) = −0.6 Với h = 0.1 ta có PHƯƠNG TRÌNH VI PHÂN Ngày 6 tháng 12 năm 2016 38 / 54 CuuDuongThanCong.com https://fb.com/tailieudientucntt Hệ phương trình vi phân Công thức Runge-Kutta bậc bốn  K1x = h.yk−1 K1y = h(−2xk−1 + 2yk−1 + e2tk−1 . sin tk−1) K2x = h ( yk−1 + K1y 2 ) K2y = hg ( tk−1 + h 2 , xk−1 + K1x 2 , yk−1 + K1y 2 ) K3x = hf ( tk−1 + h 2 , xk−1 + K2x 2 , yk−1 + K2y 2 ) K3y = hg ( tk−1 + h 2 , xk−1 + K2x 2 , yk−1 + K2y 2 ) K4x = hf (tk−1 + h, xk−1 + K3xyk−1 + K3y ), K4y = hg(tk−1 + h, xk−1 + K3x , yk−1 + K3y ) x(tk) ≈ xk = xk−1 + 1 6 (K1x + 2K2x + 2K3x + K4x) y(tk) ≈ yk = yk−1 + 1 6 (K1y + 2K2y + 2K3y + K4y ) với k = 1, 2, . . . , n. PHƯƠNG TRÌNH VI PHÂN Ngày 6 tháng 12 năm 2016 39 / 54 CuuDuongThanCong.com https://fb.com/tailieudientucntt Hệ phương trình vi phân Công thức Runge-Kutta bậc bốn tk x(tk) xk x ′(tk) yk 0.0 −0.4000 −0.4000 −0.60000 −0.6000 0.1 −0.4617 −0.4617 −0.6316 −0.6316 0.2 −0.5256 −0.5256 −0.6401 −0.6401 0.3 −0.5886 −0.5886 −0.6136 −0.6136 0.4 −0.6466 −0.6466 −0.5366 −0.5366 PHƯƠNG TRÌNH VI PHÂN Ngày 6 tháng 12 năm 2016 40 / 54 CuuDuongThanCong.com https://fb.com/tailieudientucntt Bài toán biên tuyến tính cấp 2 Đặt vấn đề Các phương pháp tìm nghiệm gần đúng của phương trình vi phân thường đòi hỏi các điều kiện được cho tại một thời điểm ban đầu nào đó. Đối với phương trình vi phân bậc hai, ta cần 2 giá trị y(x0) và y ′(x0). Tuy nhiên, nhiều bài toán trong thực tế cho thấy điều kiện của hàm cần tìm được cho tại nhiều thời điểm khác nhau. Vấn đề này dẫn tới việc tìm nghiệm gần đúng của 1 dạng bài toán thứ hai được gọi là bài toán biên. PHƯƠNG TRÌNH VI PHÂN Ngày 6 tháng 12 năm 2016 41 / 54 CuuDuongThanCong.com https://fb.com/tailieudientucntt Bài toán biên tuyến tính cấp 2 Đặt vấn đề Xét bài toán C H A P T E R 11 Boundary-Value Problems for Ordinary Differential Equations Introduction A common problem in civil engineering concerns the deflection of a beam of rectangular cross section subject to uniform loading while the ends of the beam are supported so that they undergo no deflection. x S S l w(x)0 Suppose that l, q, E, S, and I represent, respectively, the length of the beam, the intensity of the uniform load, the modulus of elasticity, the stress at the endpoints, and the central moment of inertia. The differential equation approximating the physical situation is of the form d2w dx2 (x) = S EI w(x)+ qx 2EI (x − l), where w(x) is the deflection a distance x from the left end of the beam. Since no deflection occurs at the ends of the beam, there are two boundary conditions w(0) = 0 and w(l) = 0. When the beam is of uniform thickness, the product EI is constant. In this case the exact solution is easily obtained. When the thickness is not uniform, the moment of inertia I is a function of x, and approximation techniques are required. Problems of this type are considered in Exercises 7 of Section 11.3 and 6 of Section 11.4. The differential equations in Chapter 5 are of first order and have one initial condition to satisfy. Later in the chapter we saw that the techniques could be extended to systems of equations and then to higher-order equations, but all the specified conditions are on the same endpoint. These are initial-value problems. In this chapter we show how to approximate the solution to boundary-value problems, differential equations with conditions imposed at different points. For first-order differential equations, only one condition is specified, so there is no distinction between initial-value and boundary-value problems. We will be considering second-order equations with two boundary values. Physical problems that are position-dependent rather than time-dependent are often described in terms of differential equations with conditions imposed at more than one point. 671 Copyright 2010 Cengage Learning. All Rights Reserved. May not be copied, scanned, or duplicated, in whole or in part. Due to electronic rights, some third party content may be suppressed from the eBook and/or eChapter(s). Editorial review has deemed that any suppressed content does not materially affect the overall learning experience. Cengage Learning reserves the right to remove additional content at any time if subsequent rights restrictions require it. PHƯƠNG TRÌNH VI PHÂN Ngày 6 tháng 12 năm 2016 42 / 54 CuuDuongThanCong.com https://fb.com/tailieudientucntt Bài toán biên tuyến tính cấp 2 Đặt vấn đề Trong phần này chúng ta chỉ xét bài toán biên của phương trình vi phân thường tuyến tính cấp hai với điều kiện biên được cho ở 2 điểm có dạng{ p(x)y ′′(x) + q(x)y ′(x) + r(x)y(x) = f (x), a < x < b, y(a) = α, y(b) = β. với phương pháp sai phân hữu hạn. PHƯƠNG TRÌNH VI PHÂN Ngày 6 tháng 12 năm 2016 43 / 54 CuuDuongThanCong.com https://fb.com/tailieudientucntt Bài toán biên tuyến tính cấp 2 Phương pháp sai phân hữu hạn Phương pháp sai phân hữu hạn Chọn số tự nhiên bất kỳ n > 0. Chia đều đoạn [a, b] thành n đoạn bởi các điểm chia x0 = a, xk = x0 + kh, k = 1, 2, . . . , n − 1, xn = b với h = b − a n . Tại các nút xk , k = 1, 2, . . . , n− 1 bên trong đoạn [a, b] sử dụng công thức sai phân hướng tâm, ta có y ′(xk) ≈ y(xk+1)− y(xk−1) 2h = yk+1 − yk−1 2h PHƯƠNG TRÌNH VI PHÂN Ngày 6 tháng 12 năm 2016 44 / 54 CuuDuongThanCong.com https://fb.com/tailieudientucntt Bài toán biên tuyến tính cấp 2 Phương pháp sai phân hữu hạn y ′′(xk) ≈ y(xk+1)− 2y(xk) + y(xk−1) h2 = yk+1 − 2yk + yk−1 h2 Thay vào phương trình đã cho ta được pk yk+1 − 2yk + yk−1 h2 + qk yk+1 − yk−1 2h + rkyk = fk , ∀k = 1, 2, . . . , n − 1 với pk = p(xk), qk = q(xk), rk = r(xk) và fk = f (xk). PHƯƠNG TRÌNH VI PHÂN Ngày 6 tháng 12 năm 2016 45 / 54 CuuDuongThanCong.com https://fb.com/tailieudientucntt Bài toán biên tuyến tính cấp 2 Phương pháp sai phân hữu hạn Từ các điều kiện biên y0 = α, yn = β sau khi biến đổi ta thu được hệ phương trình y0 = α, yn = β(pk h2 − qk 2h ) yk−1 + ( rk − 2pk h2 ) yk + (pk h2 + qk 2h ) yk+1 = fk ∀k = 1, 2, . . . , n − 1. Đây chính là hệ phương trình đại số tuyến tính cấp n − 1 : AY = B với A là ma trận A =  r1 − 2p1 h2 p1 h2 + q1 2h 0 . . . 0 p2 h2 − q2 2h r2 − 2p2 h2 p2 h2 + q2 2h . . . 0 . . . . . . . . . . . . . . . 0 0 0 . . . rn−1 − 2pn−1 h2  PHƯƠNG TRÌNH VI PHÂN Ngày 6 tháng 12 năm 2016 46 / 54 CuuDuongThanCong.com https://fb.com/tailieudientucntt Bài toán biên tuyến tính cấp 2 Phương pháp sai phân hữu hạn Y = [y1, y2, . . . , yn−1]T và B =  f1 − (p1 h2 − q1 2h ) α f2 . . . fn−2 fn−1 − (pn−1 h2 + qn−1 2h ) β  PHƯƠNG TRÌNH VI PHÂN Ngày 6 tháng 12 năm 2016 47 / 54 CuuDuongThanCong.com https://fb.com/tailieudientucntt Bài toán biên tuyến tính cấp 2 Ma trận 3 đường chéo Ma trận A ở trên là ma trận 3 đường chéo. Để giải hệ phương trình trên thì ta dùng phương pháp phân rã LU. A =  a11 a12 0 . . . 0 0 a21 a22 a23 . . . 0 0 0 a32 a33 . . . 0 0 . . . . . . . . . . . . . . . . . . 0 0 0 . . . an−1,n−1 an−1,n 0 0 0 . . . an,n−1 ann  PHƯƠNG TRÌNH VI PHÂN Ngày 6 tháng 12 năm 2016 48 / 54 CuuDuongThanCong.com https://fb.com/tailieudientucntt Bài toán biên tuyến tính cấp 2 Ma trận 3 đường chéo Khi đó phân rã Doolittle cho ta L =  1 0 0 . . . 0 `21 1 0 . . . 0 0 `32 1 . . . 0 . . . . . . . . . . . . . . . 0 0 0 . . . 1  , U =  u11 u12 0 . . . 0 0 u22 u23 . . . 0 0 0 u33 . . . 0 . . . . . . . . . . . . . . . 0 0 0 . . . unn  PHƯƠNG TRÌNH VI PHÂN Ngày 6 tháng 12 năm 2016 49 / 54 CuuDuongThanCong.com https://fb.com/tailieudientucntt Bài toán biên tuyến tính cấp 2 Ví dụ Ví dụ Xét bài toán biên y ′′ − y ′ − 2y = cos x , 0 < x < pi 2 y(0) = −0.3, y(pi 2 ) = −0.1 có nghiệm chính xác y(x) = −0.1(sin x + 3 cos x). Sử dụng phương pháp sai phân hữu hạn xấp xỉ nghiệm gần đúng và so sánh với nghiệm chính xác trong trường hợp h = pi 8 . PHƯƠNG TRÌNH VI PHÂN Ngày 6 tháng 12 năm 2016 50 / 54 CuuDuongThanCong.com https://fb.com/tailieudientucntt Bài toán biên tuyến tính cấp 2 Ví dụ Ta có hệ phương trình y0 = −0.3, y4 = −0.1( 1 h2 − −1 2h ) yk−1 + ( −2− 2 h2 ) yk + ( 1 h2 + −1 2h ) yk+1 = cos(xk) ∀k = 1, 2, 3. ⇔  y0 = −0.3, y4 = −0.1 ( 1 h2 + 1 2h )y0 + (−2− 2h2 )y1 + ( 1h2 − 12h)y2 = cos(pi8 ) ( 1 h2 + 1 2h )y1 + (−2− 2h2 )y2 + ( 1h2 − 12h)y3 = cos(pi4 ) ( 1 h2 + 1 2h )y2 + (−2− 2h2 )y3 + ( 1h2 − 12h)y4 = cos(3pi8 ) ⇔  (−2− 2 h2 )y1 +( 1 h2 − 1 2h )y2 +0y3 = cos( pi 8 )− ( 1 h2 + 1 2h )y0 ( 1 h2 + 1 2h )y1 +(−2− 2h2 )y2 +( 1h2 − 12h )y3 = cos(pi4 ) 0y1 +( 1 h2 + 1 2h )y2 +(−2− 2h2 )y3 = cos( 3pi8 )− ( 1h2 − 12h )y4 PHƯƠNG TRÌNH VI PHÂN Ngày 6 tháng 12 năm 2016 51 / 54 CuuDuongThanCong.com https://fb.com/tailieudientucntt Bài toán biên tuyến tính cấp 2 Ví dụ Bấm máy. pi 8 −Shift+STO+M Mode - eqn - anx+bny+cnz=dn k xk yk y(xk) |y(xk)− yk | 0 0 −0.30000 −0.30000 0.00000 1 pi 8 −0.31569 −0.31543 0.00025 2 pi 4 −0.28291 −0.28284 0.00007 3 3pi 8 −0.20700 −0.20719 0.00019 4 pi 2 −0.10000 −0.10000 0.00000 PHƯƠNG TRÌNH VI PHÂN Ngày 6 tháng 12 năm 2016 52 / 54 CuuDuongThanCong.com https://fb.com/tailieudientucntt Bài toán biên tuyến tính cấp 2 Ví dụ Bài tập 1 Sử dụng phương pháp sai phân hữu hạn giải bài toán y ′′ = 2y ′ − 2y + xeX − x 0 ≤ x ≤ 2, y(0) = 0, y(2) = −4 với h = 0.5 và so sánh sai số với nghiệm là y(x) = 16x 3ex − 53xex + 2ex − x − 2, 2 Sử dụng phương pháp sai phân hữu hạn giải bài toán y ′′ = 4(y − x), 0 ≤ x ≤ 1, y(0) = 0, y(1) = 2 với h = 1 2 và h = 1 4 và so sánh sai số với nghiệm là y(x) = e2(e4 − 1)−1(e2x − e−2x) + x , PHƯƠNG TRÌNH VI PHÂN Ngày 6 tháng 12 năm 2016 53 / 54 CuuDuongThanCong.com https://fb.com/tailieudientucntt The end THANK YOU FOR ATTENTION PHƯƠNG TRÌNH VI PHÂN Ngày 6 tháng 12 năm 2016 54 / 54 CuuDuongThanCong.com https://fb.com/tailieudientucntt