Phương pháp tính - Chương 5: Giải phương trình vi phân thường

Nội dung text: Phương pháp tính - Chương 5: Giải phương trình vi phân thường

1. BỘ MƠN TỐN ỨNG DỤNG - ĐHBK PHƯƠNG PHÁP TÍNH – SV CHƯƠNG 5 GIẢI PHƯƠNG TRÌNH VI PHÂN THƯỜNG • TS. NGUYỄN QUỐC LÂN (5/2006)
2. NỘI DUNG A- BÀI TỐN CƠSI (GIÁ TRỊ ĐẦU) 1 – PHƯƠNG PHÁP EULER 2 – EULER CẢI TIẾN + RUNGE – KUTTA 3 – HỆ PHƯƠNG TRÌNH VI PHÂN THƯỜNG 4 – PHƯƠNG TRÌNH VI PHÂN CẤP CAO B- BÀI TỐN BIÊN 1- PHƯƠNG PHÁP SAI PHÂN HỮU HẠN
3. BÀI TỐN CƠSI Tìm hàm y = y(t) thoả phương trình vi phân thường & điều kiện đầu y'= f (t, y), a t b y(a) = Giải xấp xỉ: Chia [a, b] thành n đoạn bằng nhau, độ dài h = (b – a)/n, (n + 1) điểm chia t0 = a < t1 = a + h < < tn = b y1 = ? y0 = a h b a = t0 t1 t2 b = tn Cần tính gần đúng giá trị wk yk = y(tk), k = 1 → n
4. MINH HOẠ Ý TƯỞNG y'= −5y + 5t 2 + 2t, 0 t Bài tốn Cơsi: Với bước chia h = 0.5 y(0) =1 3 f (x + h) − f (x ) & cơng thức xấp xỉ đạo hàm 2 điểm: f '(x ) 0 0 0 h hãy tính xấp xỉ nghiệm y tại t = 0.5, t = 1. Từ đĩ xây dựng đa thức nội suy Lagrange (spline) ygđ và vẽ đồ 1 thị so sánh với nghiệm chính xác g(t) = t 2 + e−5t 3 Điểm chia: t0 = 0 t1 = 0.5 t2 =1. Kết quả tìm được: 2 y(0.5) = −0.5 ygđ.Lagrange = at + bt + c 2 y(1.0) =1.875 ygđ = 6.42t − 4.87t + 0.33
5. CÁC SƠ ĐỒ GIẢI XẤP XỈ PTRÌNH VPHÂN THƯỜNG Btốn Cơsi: Tìm y(t) Sơ đồ Euler (i = 0 → n – 1) y'= f (t, y),t a,b w0 = . Giả sử wi đã biết y(a) = wi+1 = wi + hf (ti , wi ) Chia [a, b] → n đoạn S/đ Euler cải tiến (i = 0 → n – 1) b − a h = , ti = a + ih w = . Giả sử w đã biết n 0 i Tính wi, i = 0 → n k1 = hf (ti , wi ), k2 = hf (ti + h, wi + k1) wi+1 = wi + (k1 + k2 ) 2 Sơ đồ Runge – k1 = hf (ti , wi ), k2 = hf (ti + h 2, wi + k1 2) Kutta: w0 = . k3 = hf (ti + h 2, wi + k2 2), k4 = hf (ti+1, wi + k3 ) w = w + (k + 2k + 2k + k ) 6 Giả sử biết wi i+1 i 1 2 3 4
6. VÍ DỤ PHƯƠNG PHÁP EULER Bằng p/pháp Euler, giải bài tốn Cơsi với n = 3 đoạn chia: y'= y − t2 +1, 0 t 1 y(0) = 0.5 So sánh nghiệm xấp xỉ với nghiệm g(t) = (t+1)2 – 0.5et. Từ 1 đĩ tính xấp xỉ tích phân bằng c/t hình thang: I = y(t)dt 0 Giải: f(t,y) = y – t2 + 1 t0 = 0, w0 = 0.5 t1, w1 t2 , w2 t3, w3 h = (b–a)/n = 1/3 w = 0.5 Sơ đồ Euler: 0 2 wi+1 = wi + hf (ti , wi ) = wi + 0.2(wi −ti +1)
7. KẾT QUẢ PHƯƠNG PHÁP EULER Bảng kết quả: i ti wi gi = g(ti ) | gi - wi | 0 0 0.5 0.5 0 1 1/3 2 2/3 3 1. Tính gần đúng tích phân với cơng thức hình thang 1 h h y(t)dt y t + 2y t + 2y t + y t w + 2w + 2w + w  ( 0 ) ( 1 ) ( 2 ) ( 3 )  0 1 2 3  0 2 2 =1.3528807
8. VÍ DỤ EULER CẢI TIẾN Tính y(1.) của bt Cơsi sau bằng y'= y − t2 +1, 0 t 1 SĐ Euler cải tiến với h = 0.5: y(0) = 0.5 f (t, y) = y −t 2 +1, h = 0.5 t0 = 0, = 0.5 t1 = 0.5, w1 = ? t2 , w2 k + k k = hf (t , w ) , k = hf (t + h, w + k ) → w = w + 1 2 1 i i 2 i i 1 i+1 i 2 i ti wi k1 k2 0 0.0 0.5 0.75 1.0 1 0.5 1.375 1.0625 1.21875 2 1.0 2.515625
9. VÍ DỤ RUNGE – KUTTA y'= y − t 2 +1 Tính y(1.) bằng Runge – Kutta với h = 0.5 y(0) = 0.5 h k k = hf (t , w ) , k = hf (t + , w + 1 ) Runge – 1 i i 2 i 2 i 2 Kutta 4: k3 = hf (ti + h 2, wi + k2 2) , k4 = hf (ti + h, wi + k3 ) w → w k + 2k + 2k + k i i+1 w = w + 1 2 3 4 i+1 i 6 i ti wi k1 k2 k3 k4 0 0.0 0.5 0.75 0.90625 0.9451325 1.0976563 1 0.5 1.4251302 1.0875651 1.2032064 1.2331167 1.3286235 2 1.0 2.6396027
10. HỆ PHƯƠNG TRÌNH VI PHÂN THƯỜNG Bài tốn Cơsi : Tìm hai hàm u1 = u1(t), u2 = u2(t) thoả du1 = f1(t,u1,u2 ), a t b dt u1(a) = 1 & Điều kiện đầu du u (a) = 2 = f (t,u ,u ), a t b 2 2 dt 2 1 2 Chia [a, b] thành đoạn bằng nhau: Phân hoạch & rời rạc hố 0 0 1 1 w1 = 1, w2 = 2 w1 , w2 1, 2 2 2 u1(t1 ), u2 (t1 ) = ? w1 , w2 = ? a = t0 t1 = a + h t2 = a + 2h i i 0 0 i i Ký hiệu: w1 u1(ti ), w2 u2 (ti ), i 0 Biết w1 , w2 tính w1,w2 ?
11. MINH HOẠ Ý TƯỞNG Xét bài tốn Cơsi với hệ phương trình vi phân thường: 2 2t u1'= 3u1 + 2u2 − (2t +1)e , u1(0) =1 2 2t u2 '= 4u1 + u2 + (t + 2t − 4)e , u2 (0) =1 Với bước chia h = 0.5, tính xấp xỉ nghiệm u1, u2 tại t = 0.5; 1 So sánh giá trị tính được với giá trị nghiệm chính xác: 1 1 1 2 u (t) = e5t − e−t + e2t ; u (t) = e5t + e−t + t 2e2t 1 3 3 2 3 3 t0 = 0 t1 = 0.5 t2 =1. Điểm chia: t u1 u2 0 1 1 0.5 u1(0) =1 u1(0.5) = ? u1(1) = ? 1.0 u2 (0) =1 u2 (0.5) = ? u2 (1) = ?
12. SƠ ĐỒ EULER Bài tốn Cơsi : Tìm hai hàm u1 = u1(t), u2 = u2(t) thoả du1 = f1(t,u1,u2 ), a t b dt u1(a) = 1 & Điều kiện đầu du u (a) = 2 = f (t,u ,u ), a t b 2 2 dt 2 1 2 0 0 i i S/đồ Euler: w1 = 1, w2 = 2 Giả sử biết w1 , w2 (i = 0 → n −1) i+1 i i i i+1 i i i w1 = w1 + hf1(ti , w1 , w2 ), w2 = w2 + hf2 (ti , w1 , w2 ) 0 2 2t =1 w =1 u '= 3u + 2u − 2t +1 e , u (0) =1 1 1 1 12 () 1 =1 0 f1(t,u1,u2 ) 2 w2 =1 VD: u '= 4u + u + t 2 + 2t − 4 e2t , u (0) =1 w1 =1+ 0.5 f (0,1,1) = 2 12 () 2 1 1 1 f2 (t,u1,u2 ) w2 =1+ 0.5 f2 (0,1,1) =
13. ÁP DỤNG : PHƯƠNG TRÌNH VI PHÂN CẤP 2 Bài tốn Cơsi cấp 2 (Ph/trình vi phân cấp 2 và đkiện đầu): y"= f (t, y, y'), t a y(a) = 1, y'(a) = 2 Đưa về bài tốn Cơsi cấp 1: Đổi biến u1(t)= y(t), u2(t)=y’(t) u1'= u2 = f1(t,u1,u2 ) u2 '= y''= f (t, y, y') = f2 (t,u1,u2 ) 0 u1(a) = y(a) = 1 w1 = 1 Điều kiện đầu: u (a) = y'(a) = 0 2 2 w2 = 2 w1 = w0 + hf (t , w0 , w0 )= + h Sơ đồ Euler: 1 1 1 0 1 2 1 2 1 0 0 0 w2 = w2 + hf2 (t0 , w1 , w2 )= 2 + hf2 (a, 1, 2 )
14. VÍ DỤ Với h = 0.1, tính xấp xỉ giá trị y(0.2), y’(0.2) của nghiệm bài tốn sau bằng phương pháp Euler: y"−2y'+2y = e2t sin t , t 0 y(0) = −0.4, y'(0) = −0.6 Đổi biến đưa về bài tốn Cơsi cấp 1: u1 = y(t), u2 = y’(t) u1'= u2 = f1(t,u1,u2 ) & u1(0) = −0.4 , u2 (0) = −0.6 2t 2t u2 '= −2y + 2y'+e sin t = −2u1 + 2u2 + e sin t = f2 (t,u1,u2 ) w0 = −0.4 w1 = w0 + hf (t , w0 , w0 )= −0.4 + 0.1f (0,−0.4,−0.6) 1 1 1 1 0 1 2 1 0 1 0 0 0 w2 = −0.6 w2 = w2 + hf2 (t0 , w1 , w2 )= −0.6 + 0.1f2 (0,−0.4,−0.6)
15. BÀI TỐN BIÊN Bài tốn biên cấp 2: Tìm hàm y = y(x) thoả phương trình y''= f (x, y, y'), a x b y(a) = , y(b) =  Hay gặp: Bài tốn biên tuyến tính cấp 2 y"= p(x)y'+q(x)y + r(x), a x b y(a) = , y(b) = 
16. MINH HOẠ Tính giá trị nghiệm y của bài tốn biên tuyến tính cấp 2 y''= −(x +1)y'+2y + (1− x2 )e−x , 0 x 1 y(0) = −1, y(1) = 0 tại các điểm chia cách đều của [0, 1] với bước chia h = 1/3 và xấp xỉ đạo hàm y’, y’’ bằng cơng thức hướng tâm Điểm chia: x0 = 0 x1 =1 3 x2 = 2 3 x3 =1 1 2 y(0) = y = −1 y = y = ? y = y = ? y(1) = y = 0 0 3 1 3 2 3
17. PHƯƠNG PHÁP SAI PHÂN HỮU HẠN y"= p(x)y'+q(x)y + r(x), a x b (*) BT biên tuyến tính y(a) = , y(b) =  Chia [a, b] thành các đoạn nhỏ bằng nhau. Thay x = xk vào (*). Xấp xỉ y’(xk) , y’’(xk): cơng thức đạo hàm hướng tâm y(x )− y(x ) y(x )− y(x ) y'(x ) 2 0 y'(x ) 3 1 1 2h 2 2h h a= x0 x1 x2 x3 b= xn+1 y(x )− 2y(x )+ y(x ) y(x )− 2y(x )+ y(x ) y"(x ) 2 1 0 y"(x ) 3 2 1 1 h2 2 h2
18. CƠNG THỨC LẮP GHÉP n mốc xk (a, b) – ứng n giá trị yk chưa biết → Ma trận cấp n T Ký hiệu pk = p(xk) yk = y(xk), 1 k n y= [y1, yn] : Ay = b 2 h 2 h 2 + h q1 −1+ p1 0  0 − h r1 + 1+ p1 2 2 h 2 h 2 −1− p 2 + h q −1+ p   − h r2 2 2 2 2 2 − h2r A = 0    0 b = 3  h     −1+ pn−1 2 2 − h rn−1 h 2 2 h 0  0 −1− pn 2 + h qn − h rn + 1− pn  2 2
19. LẬP BẢNG LẮP GHÉP y"= p(x)y'+q(x)y + r(x), a x b (*) BT biên tuyến tính y(a) = , y(b) =  ❖ Chia [a, b] thành các đoạn nhỏ độ dài h. n điểm chia xk (khơng kể 2 đầu) – ứng với yk chưa biết → n ẩn số yk ❖ Lập bảng cột xk → pk = p(xk), qk = q(xk), rk = r(xk) → akk (đ/chéo chính), ak,k+1 (chéo trên), ak-1,k (dưới), bk → Nghiệm yk ❖ Đ/chéo akk: k = 1 → n; ak,k+1: k = 1 → (n – 1), ak-1,k: k = 2 → n i xk pk qk rk akk ak,k+1 ak-1,k bk yk 1 2 3
20. VÍ DỤ Giải bài tốn biên cấp 2 sau bằng phương pháp sai phân hữu hạn với bước chia h = 0.2 y"= −3y'+2y + 2x + 3 y(0) = 2, y(1) =1 h = 0.2 n = 5 6 điểm chia Hệ phương trình 4 ẩn Ma trận cấp 4: Chéo chính akk – 4 phần tử; Chéo trên ak, k+1: 3 i xi pi qi ri aii ai,i+1 ai-1,i bi 1 0.2 − 3 2 3.4 − 2.08 1.3 −1.264 2 0.4 3.8 1.3 0.7 0.152 3 0.6 4.2 0.7 0.168 4 0.8 4.6 −1.116
21. KẾT QUẢ Giải hệ bằng phép khử Gauss (làm trịn 3 chữ số lẻ): 2.08 −1.3 0 0 1.264 − 0.7 2.08 −1.3 0 − 0.152 Ab= 0 − 0.7 2.08 −1.3 − 0.168 0 0 − 0.7 2.08 1.116 1 − 0.625 0 0 0.608 1.006 0 1.642 −1.3 0 0.273 0.636 → Ab= y = 0 − 0.7 2.08 −1.3 − 0.168 0.593 0 0 − 0.7 2.08 1.116 0.736