Bài giảng Phương pháp số - Bài 6: Nghiệm của các phương trình vi phân

Nội dung text: Bài giảng Phương pháp số - Bài 6: Nghiệm của các phương trình vi phân

1. BÀI 6 NGHIỆM CỦA CÁC PHƢƠNG TRÌNH VI PHÂN
2. MỘT SỐ KHÁI NIỆM CƠ SỞ (1) 1. PHƢƠNG TRÌNH VI PHÂN (PTVP) CẤP n y(n)(x) = f(x, y(x), y´(x) , , y(n–1)(x)) 2. NGHIỆM CỦA PTVP CẤP n: Nghiệm của PTVP cấp n là một hàm Φ(x) khả vi liên tục đến cấp n trên khoảng đang xét và Φ(n)(x) = f(x, Φ(x), Φ’(x), Φ’’(x), , Φ(n-1)(x)) Nghiệm tổng quát của PTVP cấp n thƣờng chứa n hằng số tùy ý, và do vậy tồn tại một họ n tham số của các nghiệm Nghiệm của bài toán giá trị đầu: Tìm nghiệm riêng suy từ nghiệm tổng quát của PTVP sao cho nó thỏa mãn n giá trị đầu đã biết (n–1) y(x0), y’(x0), , y (x0) PHƢƠNG PHÁP SỐ-Bài 6 2
3. MỘT SỐ KHÁI NIỆM CƠ SỞ (2) 3. NGHIỆM CỦA PTVP TUYẾN TÍNH CẤP n: • PTVP là tuyến tính cấp n (n) (n–1) y (x) + p1(x)y (x) + + pn-1(x)y´(x) + pn(x)y(x) = q(x) • PTVP tuyến tính thuần nhất cấp n (khi vế phải q(x) = 0) i> Nếu y1(x), y2(x), . . . , ym(x) là các nghiệm bất kì của PTVP tuyến tính thuần nhất thì tổ hợp tuyến tính của chúng C1y1(x) + C2y2(x) + · · · + Cmym(x) cũng là một nghiệm Ii> Các nghiệm y1(x), y2(x), . . . , ym(x) này đƣợc gọi là độc lập tuyến tính nếu định thức hàm ' (m 1) y1 y1 y1 y y' y(m 1) 2 2 2 0 ' (m 1) ym ym ym PHƢƠNG PHÁP SỐ-Bài 6 3
4. MỘT SỐ KHÁI NIỆM CƠ SỞ (3) 3. NGHIỆM CỦA PTVP TUYẾN TÍNH CẤP n: • PTVP tuyến tính thuần nhất cấp n (tiếp) iii> Nếu y1(x), y2(x), . . . , yn(x) là n nghiệm độc lập tuyến tính của PTVP thuần nhất cấp n, thì tổ hợp tuyến tính của chúng Y(x) = C1y1(x) + C2y2(x) + · · · + Cnyn(x) đƣợc gọi là nghiệm tổng quát. • PTVP tuyến tính thuần nhất cấp n với các hệ số hằng (n) (n–1) y + an–1y + · · · + a0y = 0 (*) có nghiệm dạng eβx với β thỏa mãn PT đặc tính n n–1 β + an–1β + · · · + a0 = 0 ( ) PHƢƠNG PHÁP SỐ-Bài 6 4
5. MỘT SỐ KHÁI NIỆM CƠ SỞ (4) 3. NGHIỆM CỦA PTVP TUYẾN TÍNH CẤP n: • PTVP tuyến tính thuần nhất cấp n với các hệ số hằng: - Nếu phƣơng trình ( ) có n nghiệm thực phân biệt βi thì β1x β2x βnx Y(x) C1e C2e Cne là nghiệm tổng quát của của PTVP thuần nhất cấp n (*) - Nếu β1 = α + iβ là một nghiệm phức của ( ), thì β2 = α – iβ cũng là một nghiệm của PT đặc tính ( ) - Nếu β là một nghiệm kép của PT đặc tính ( ), thì βx βx y1 e và y2 xe cũng là các nghiệm độc lập tuyến tính của PT thuần nhất (*) PHƢƠNG PHÁP SỐ-Bài 6 5
6. MỘT SỐ KHÁI NIỆM CƠ SỞ (5) 3. NGHIỆM CỦA PTVP TUYẾN TÍNH CẤP n: • PTVP tuyến tính không thuần nhất cấp n hệ số hằng: (n) (n–1) y + an–1y + · · · + a0y = q(x) ( ) Nếu PTVP ( ) có một nghiệm riêng ζ(x), thì nghiệm tổng quát của PTVP tuyến tính không thuần nhất là y(x) ζ(x) Y(x) β1x β2x βnx ζ(x) C1e C2e Cne với giả thiết rằng βi là các nghiệm thực phân biệt của PTVP đặc tính ( ) tƣơng ứng PHƢƠNG PHÁP SỐ-Bài 6 6
7. MỘT SỐ KHÁI NIỆM CƠ SỞ (6) 3. NGHIỆM CỦA PTVP TUYẾN TÍNH CẤP n: • PTVP tuyến tính không thuần nhất cấp n hệ số hằng: Ví dụ: Tìm nghiệm tổng quát của phƣơng trình sau y´´ – 4y´ + 3y = x 2 Giải: Từ PT đặc tính β - 4β + 3 = 0 β1 = 1, β2 = 3 nghiệm tổng quát của PT thuần nhất là x 3x Y(x) C1e C2e Ta tìm một nghiệm riêng y1(x) = ax + b của PT ban đầu a = 1/3, b = 4/9 Vậy nghiệm tổng quát của PT không thuần nhất là 1 4 y(x) y (x) Y(x) x C ex C e3x 1 3 9 1 2 PHƢƠNG PHÁP SỐ-Bài 6 7
8. MỘT SỐ KHÁI NIỆM CƠ SỞ (7) 4. BÀI TOÁN CÔ SI CỦA PTVP CẤP n Tìm nghiệm riêng của PTVP (phi tuyến) cấp n y(n) = f(x, y, y’, , y(n-1)) thỏa mãn các giá trị đầu (n 1) (n 1) y(x0 ) y0 , y'(x0 ) y'0 , , y (x0 ) y0 Ví dụ 1: Tìm nghiệm riêng của PTVP y’ = xy thỏa mãn đ/k giá trị đầu y(0) = 1 Giải: dy . dy xy x dx tích phân 2 vế ta có nghiệm TQ dx y dy x2 x2 x dx ln | y | C y Ae 2 y 2 x2 Từ đ/k y(0) = 1 ta tìm đƣợc nghiệm riêng y e 2 PHƢƠNG PHÁP SỐ-Bài 6 8
9. MỘT SỐ KHÁI NIỆM CƠ SỞ (8) 5. THUẬT TOÁN TAYLOR Giải BT Côsi: Tìm nghiệm của PTVP cấp một y´ = f(x,y) biết giá trị đầu y(x0) = y0 ( ) • Giả sử f là đủ khả vi theo cả x và y, ( ) có nghiệm duy nhất nếu ∂f/∂y liên tục trên miền đang xét. Khai triển nghiệm y(x) thành chuỗi Taylor xung quanh điểm x0: (x x )2 y(x) y (x x )y'(x ) 0 y''(x )  0 0 0 2! 0 • Tính các đạo hàm bằng cách lấy đạo hàm toàn phần của ( ) theo x: y´ = f(x,y) y´´ = f´ = fx + fyy´ = fx + fyf 2 2 y´´´ = f´´ = fxx + fxyf + fyxf + fyyf + fyfx + fy f 2 2 = fxx + 2fxyf + fyyf + fxfy + fy f PHƢƠNG PHÁP SỐ-Bài 6 9
10. MỘT SỐ KHÁI NIỆM CƠ SỞ (9) 5. THUẬT TOÁN TAYLOR (tiếp) Thuật toán Taylor bậc k tìm một nghiệm xấp xỉ của BT Côsi y´ = f(x,y) thỏa mãn đ/k đầu y(a) = y0 trên [a, b] 1. Chọn bƣớc h = (b – a) / N, đặt xn = a + nh, n = 0, 1, . . . , N 2. Tạo ra các xấp xỉ yn cho y(xn) từ phép đệ quy yn+1 = yn + hTk(xn,yn ) n = 0, l, . . . , N – 1 trong đó Tk(x, y) đƣợc định nghĩa bởi h h k 1 T (x, y) f (x, y) f'(x, y)  f (k 1)(x, y), k 1,2, k 2! k! Sai số địa phƣơng của thuật toán Taylor bậc k là h k 1f (k) (, y( )) h k 1 E y(k 1) ( ) , x  x h (k 1)! (k 1)! n n PHƢƠNG PHÁP SỐ-Bài 6 10
11. MỘT SỐ KHÁI NIỆM CƠ SỞ (10) 5. THUẬT TOÁN TAYLOR (tiếp) Ví dụ: tìm một nghiệm xấp xỉ của BT Côsi 1 y y' y2 tháa mãn đ / k y(1) 1 x 2 x trên [1, 2] bằng thuật toán Taylor bậc 2 với h = 0.1 Đặt xi = a + ih, i = 0, 1, . . . , 10 1 y 2 y' y f (x, y) y2 , f'(x, y) 2yy' x 2 x x 3 x x 2 h do đó T (x, y) f f' 2 2 h và yi 1 yi h f (xi , yi ) f'(xi , yi ) 2 thuật toán Taylor bậc cao hơn khó đƣợc chấp nhận PHƢƠNG PHÁP SỐ-Bài 6 11
12. MỘT SỐ KHÁI NIỆM CƠ SỞ (11) 5. THUẬT TOÁN TAYLOR (tiếp) Ví dụ: (tiếp) i xi yi y'=f(xi,yi) hf(xi,yi) f'(xi,yi) h^2*f'(x,y)/2 0 1 -1 1 0.1 -2 -0.010000 1 1.1 -0.910000 0.825619 0.082562 -1.502632 -0.007513 2 1.2 -0.834951 0.693094 0.069309 -1.157414 -0.005787 3 1.3 -0.771429 0.590020 0.059002 -0.910343 -0.004552 4 1.4 -0.716979 0.508273 0.050827 -0.728879 -0.003644 5 1.5 -0.669796 0.442349 0.044235 -0.592612 -0.002963 6 1.6 -0.628524 0.388410 0.038841 -0.488305 -0.002442 7 1.7 -0.592124 0.343718 0.034372 -0.407110 -0.002036 8 1.8 -0.559788 0.306273 0.030627 -0.342966 -0.001715 9 1.9 -0.530876 0.274588 0.027459 -0.291621 -0.001458 10 2 -0.504875 0.247539 0.024754 -0.250036 -0.001250 PHƢƠNG PHÁP SỐ-Bài 6 12
13. CÁC PHƢƠNG PHÁP SỐ GIẢI BT CÔSI (1) 1. PHƢƠNG PHÁP EULER GIẢI BT CÔSI CẤP 1 • Áp dụng thuật toán Taylor với k = 1 trên [a, b] yi 1 yi hf(xi , yi ) h 2 có sai số địa phƣơng E y''( ) 2 sai số toàn cục trên [a, b] cùng bậc với h cần h rất nhỏ • Ví dụ, giải BT Côsi y´ = y biết giá trị đầu y(0) = 1, h = 0.01 i xi yi f(x,y) hf(x,y) 0 0 1 1.000000 0.010000 1 0.01 1.010000 1.010000 0.010100 2 0.02 1.020100 1.020100 0.010201 3 0.03 1.030301 1.030301 0.010303 4 0.04 1.040604 1.040604 0.010406 5 0.05 1.051010 1.051010 0.010510 PHƢƠNG PHÁP SỐ-Bài 6 13
14. CÁC PHƢƠNG PHÁP SỐ GIẢI BT CÔSI (2) void Euler (double a, double b, double h, double y0) { int n = ceil((b-a)/h); vector x(n+1,0), y(n+1,0); x[0] = a; y[0] = y0; for (int i = 1; i <= n; i++) { x[i] = x[i-1] + h; y[i] = y[i-1] + h * f(x[i-1], y[i-1]); } cout << "Nghiem gan dung tren doan [" << x[0] << "," << x[n] << "]" << endl; cout << "\t i \t x \t y \t f(x,y) \t hf(x,y)" << endl; for (int i = 0; i <= n; i++) { cout << "\t "<< i << "\t " << x[i] << "\t " << y[i] << "\t " << f(x[i], y[i]) << "\t " << h*f(x[i], y[i]) << endl; } } PHƢƠNG PHÁP SỐ-Bài 6 14
15. CÁC PHƢƠNG PHÁP SỐ GIẢI BT CÔSI (3) 2. PHƢƠNG PHÁP RUNGE-KUTTA BẬC 2 • Đặt vấn đề: Đối với nghiệm của BT Côsi của PTVP cấp 1 Sử dụng CT y’ = f(x, y), y(x0) = y0 Euler, tính độ Ƣớc lƣợng hàm f(x, y) tại các điểm nghiêng ở xi+1 đƣợc lựa chọn trên mỗi khoảng con. yi+l = yi + ak1(h) + bk2(h) trong đó a = b = 1/2 k1 = hf(xi, yi) k2 = hf(xi + h, yi + k1) yi+l = yi + 1/2(k1 + k2) • Sai số địa phƣơng h3 y(x ) y (f 2ff f f 2 2f f f 2f ) O(h 4 ) i 1 i 1 12 xx xy yy x y y PHƢƠNG PHÁP SỐ-Bài 6 15
16. CÁC PHƢƠNG PHÁP SỐ GIẢI BT CÔSI (4) 2. PHƢƠNG PHÁP RUNGE-KUTTA BẬC 2 • Ví dụ: giải PTVP y'=1-y/x biết y(2) = 2, h=0.01trên [2; 2.1] i xi yi f(xi,yi) hf(xi,yi) zi f(xi+1,zi) hf(xi+1,zi) 0 2.00 2.000000 0.000000 0.000000 2.000000 0.004975 4.9751E-05 1 2.01 2.000025 0.004963 0.000050 2.000075 0.009864 9.8641E-05 2 2.02 2.000099 0.009852 0.000099 2.000198 0.014681 0.00014681 3 2.03 2.000222 0.014669 0.000147 2.000368 0.019427 0.00019427 4 2.04 2.000392 0.019416 0.000194 2.000586 0.024104 0.00024104 5 2.05 2.000610 0.024093 0.000241 2.000851 0.028713 0.00028713 6 2.06 2.000874 0.028702 0.000287 2.001161 0.033256 0.00033256 7 2.07 2.001184 0.033245 0.000332 2.001516 0.037733 0.00037733 8 2.08 2.001538 0.037722 0.000377 2.001916 0.042146 0.00042146 9 2.09 2.001938 0.042135 0.000421 2.002359 0.046496 0.00046496 10 2.10 2.002381 0.046485 0.000465 2.002846 PHƢƠNG PHÁP SỐ-Bài 6 16
17. CÁC PHƢƠNG PHÁP SỐ GIẢI BT CÔSI (5) void Runge_Kutta2 (double h, double x, double y, unsigned n) { double xs, z; cout << "Nghiem cua phuong trinh " << endl; cout << "i x y f(x,y) f(x,z)"; for (unsigned i = 0; i <= n; i++) { z = y + h * f(x, y); xs = x + h; cout << "\n“ << i << "\t “ << xs << "\t “ << y << "\t “ << f(x, y) << "\t “ << f(xs, z) << endl; y = y + h * (f(x, y) + f(xs, z)) / 2; x = xs; } } PHƢƠNG PHÁP SỐ-Bài 6 17
18. CÁC PHƢƠNG PHÁP SỐ GIẢI BT CÔSI (6) 3. PHƢƠNG PHÁP RUNGE-KUTTA BẬC 4 Đối với PTVP y’ = f(x, y), y(x0) = y0, tạo xấp xỉ cho y(x0+h) y(x0+h) = y(x0) + r1k1(h) + r2k2(h) + r3k3(h) + r4k4(h) + φ4(h) hay y(x0+h) = y0 + ∆y0 + φ4(h) trong đó ∆y0 = r1k1(h) + r2k2(h) + r3k3(h) + r4k4(h) ki(h) = hf(ξi, ૪i); ξi = x0 + ai h, a1 = 0 ૪i = y0 + bi1k1(h) + + bi, i-1 ki-1(h), i = 1, 2, 3, 4 φ4(h) = y(x0+h) - y0 - ∆y0 là sai số (địa phƣơng ) Các hệ số ai, bij, ri đƣợc chọn sao cho các đạo hàm cấp i của φ4(h) thoả mãn điều kiện (i) (i 1) 4 (0) 0 (i 0,4) và 4 (0) 0i 4 PHƢƠNG PHÁP SỐ-Bài 6 18
19. CÁC PHƢƠNG PHÁP SỐ GIẢI BT CÔSI (7) 3. PHƢƠNG PHÁP RUNGE-KUTTA BẬC 4 hệ 11 phƣơng trình, 13 ẩn a2, a3, a4, b21, b31, b32, b41, b42, b43, r1, r2, r3, r4. • Công thức Runge – Kutta bậc 4 : 1 y y k 2k 2k k i 1 i 6 1 2 3 4 k hf(x , y ) trong đó 1 i i 1 1 k 2 hf x i h, yi k1 2 2 1 1 k 3 hf x i h, yi k 2 2 2 k 4 hf(xi h, yi k 3) PHƢƠNG PHÁP SỐ-Bài 6 19
20. CÁC PHƢƠNG PHÁP SỐ GIẢI BT CÔSI (8) 3. PHƢƠNG PHÁP RUNGE-KUTTA BẬC 4 • Sai số địa phƣơng của công thức là O(h4) • Ƣu điểm: Tự bắt đầu (chỉ đòi hỏi giá trị y(x0) = y0) • Nhƣợc điểm: Mỗi bƣớc đòi hỏi 4 ƣớc lƣợng hàm y k2 yi+1 k k 3 y 1 i k4 ktổ hợp x x x i xi+1/2 i+1 PHƢƠNG PHÁP SỐ-Bài 6 20
21. CÁC PHƢƠNG PHÁP SỐ GIẢI BT CÔSI (9) 3. PHƢƠNG PHÁP RUNGE-KUTTA BẬC 4 Ví dụ: Tìm nghiệm của BT Côsi y´ = x – 1/y với điều kiện đầu y(1) = 2 trên [1; 1.4], h = 0.2 i xi yi k1=hf(xi,yi) k2=hf(xi+h/2,yi+k1/2) k3=hf(xi+h/2,yi+k2/2) k4=hf(xi+h,yi+k3) Delta 0 1 2 0.1 1.1 2.050000 0.122439 1.1 2.061220 0.122970 1.2 2.122970 0.145792 0.122768 y1 = y(1.2) ≈ y(1) + Δy0 = 2.122768 1 1.2 2.122768 0.145783 1.3 2.195660 0.168911 1.3 2.207224 0.169388 1.4 2.292157 0.192746 0.169188 y2 = y(1.4) ≈ y(1.2) + Δy1 = 2.291957 PHƢƠNG PHÁP SỐ-Bài 6 21
22. CÁC PHƢƠNG PHÁP SỐ GIẢI BT CÔSI (10) void Runge_Kutta4 (double h, double x, double y, unsigned n) { double k1, k2, k3, k4, delta; cout << " Nghiem cua phuong trinh y'= 1+y/x" << endl; cout << "\t i x y delta“ << endl; for (unsigned i = 1; i <= n; i++) { k1 = h * f(x, y); k2 = h * f(x + h / 2, y + k1 / 2); k3 = h * f(x + h / 2, y + k2 / 2); k4 = h * f(x+h, y+k3); delta= (k1 + 2 * k2 + 2 * k3 + k4) / 6; x += h; y += delta; cout << "\t " << i << "\t " << x << "\t " << y << "\t " << delta << endl; } } PHƢƠNG PHÁP SỐ-Bài 6 22
23. CÁC PHƢƠNG PHÁP SỐ GIẢI BT CÔSI (11) 4. BÀI TOÁN CÔSI CỦA HỆ PTVP CẤP 1 Các phƣơng pháp trên có thể áp dụng đƣợc cho hệ PTVP cấp một. Chẳng hạn giải bài toán Côsi sau đây: Cho hệ hai PTVP y' f (x, y,z) z' g(x, y,z) với điều kiện đầu y(x0) = y0 và z(x0) = z0, bằng phƣơng pháp Runge – Kutta bậc 4. Khi đó công thức Runge – Kutta bậc 4 có dạng yi+1≈yi+(k1+2k2+2k3+k4) / 6 zi+1≈zi+(l1+2l2+2l3+l4) / 6 PHƢƠNG PHÁP SỐ-Bài 6 23
24. CÁC PHƢƠNG PHÁP SỐ GIẢI BT CÔSI (12) 4. BÀI TOÁN CÔSI CỦA HỆ PTVP CẤP 1 Các hệ số ki và li đƣợc tính nhƣ sau k1 = hf(xi, yi, zi) l1 = hg(xi, yi, zi) k2= hf(xi+h/2, yi+k1/2, zi+l1/2) l2= hg(xi+h/2, yi+k1/2, zi+l1/2) k3 = hf(xi+h/2, yi+k2/2, zi+l2/2) l3 = hg(xi+h/2, yi+k2/2, zi+l2/2) k4 = hf(xi+1, yi+k3, zi+l3) l4 = hg(xi+1, yi+k3, zi+l3) Ví dụ: Giải BT Côsi của hệ dƣới đây trên [0; 0,4] với h = 0,1 y' z z' x zx y với các điều kiện đầu y(0) = 0 và z(0) = 1 PHƢƠNG PHÁP SỐ-Bài 6 24
25. CÁC PHƢƠNG PHÁP SỐ GIẢI BT CÔSI (13) 4. BÀI TOÁN CÔSI CỦA HỆ PTVP CẤP 1 i x y z k=hz l = h(x+xz−y) 0 0 0 1 0,1 0 0,05 0,05 1 0,1 0.005 0,05 0,05 0,0025 0,10025 0.005013 0,1 0.10025 1.005013 0,100501 0.010025 y1 = y0 + ∆y0 = 0 + 0,100167 = 0.100167 z1 = z0 + ∆z0 = 1 + 0.005008 = 1.005008 1 0,1 0,100167 1.005008 0,100501 0.010033 0,15 0,150167 1.010025 0,101003 0.015134 0,15 0,150668 1.012575 0,101258 0.015122 0,20 0,201424 1.020130 0,102013 0.020260 y2 = y1 + ∆y1 = 0.100167 + 0.101172 = 0.201339 z2 = z1 + ∆z1 = 1.005008 + 0.015134 = 1.020142 PHƢƠNG PHÁP SỐ-Bài 6 25
26. CÁC PHƢƠNG PHÁP SỐ GIẢI BT CÔSI (14) 5. BÀI TOÁN CÔSI CỦA PTVP CẤP CAO Đƣa về bài toán Cô si của hệ PTVP cấp một Ví dụ: Giải BT Côsi của phƣơng trình cấp 2 y´´ – y´x + y = x thỏa mãn các điều kiện đầu y(0) = 0 và y´(0) = 1 Cách giải: Đƣa về BT Côsi của hệ 2 PTVP cấp một y' z tƣơng đƣơng z' x zx y với các đ/k đầu y(0) = 0 và z(0) = 1. Giải hệ này tìm được y PHƢƠNG PHÁP SỐ-Bài 6 26
27. BÀI TẬP 1. Bài tập tính toán: 8.1-1 (a, b, d), 8.1-2, 8.3-3, 8.4-1 (b và d), 8.5-1, 8.12-1. Dùng phƣơng pháp Runge –Kutta bậc 4 giải lại bài tập 8.1-2 (a, c) trên [0; 0.3] với h = 0.1 2. Bài tập lập trình: Phƣơng pháp Runge –Kutta bậc 2 và bậc 4 giải BT Côsi của hệ 2 PTVP cấp 1 PHƢƠNG PHÁP SỐ-Bài 6 27