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 - Nguyễn Thị Vinh

BÀI 6 NGHIỆM CỦA CÁC PHƢƠNG TRÌNH VI PHÂN 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 y(x0), y’(x0), , y (n–1)(x0) PHƢƠNG PHÁP SỐ-Bài 6 2 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 y(n)(x) + p1(x)y (n–1)(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 0 yyy yyy yyy 1m mmm 1m 222 1m 111     )(' )(' )(' ... ............ ... ... PHƢƠNG PHÁP SỐ-Bài 6 3 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 y(n) + an–1y (n–1) + · · · + a0y = 0 (*) có nghiệm dạng eβx với β thỏa mãn PT đặc tính βn + an–1β n–1 + · · · + a0 = 0 (**) PHƢƠNG PHÁP SỐ-Bài 6 4 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ì 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ì cũng là các nghiệm độc lập tuyến tính của PT thuần nhất (*) xnβ n xβ 2 xβ 1 eC...eCeCY(x) 21  βx 2 βx 1 xeyvàey  PHƢƠNG PHÁP SỐ-Bài 6 5 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: y(n) + an–1y (n–1) + · · · + 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à 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 xnβ n xβ 2 xβ 1 eC...eCeCζ(x) xYζ(x)xy 21   )()( PHƢƠNG PHÁP SỐ-Bài 6 6 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 Giải: Từ PT đặc tính β2 - 4β + 3 = 0  β1 = 1, β2 = 3  nghiệm tổng quát của PT thuần nhất là 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à x3 2 x 1 eCeCxY )( PHƢƠNG PHÁP SỐ-Bài 6 7 x3 2 x 11 eCeC 9 4 x 3 1 xYxyxy  )()()( 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 PHƢƠNG PHÁP SỐ-Bài 6 8 )()( )(,,')(',)( 1n 00 1n 0000 yxyyxyyxy    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: . dxx y dy xy dx dy  tích phân 2 vế ta có nghiệm TQ 2 2x AeyC 2 x |y|lndxx y dy 2   Từ đ/k y(0) = 1 ta tìm đƣợc nghiệm riêng 2 2x ey  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: • 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 y´´´ = f´´ = fxx + fxyf + fyxf + fyyf 2 + fyfx + fy 2f = fxx + 2fxyf + fyyf 2 + fxfy + fy 2f    )('' ! )( )(')()( 0 2 0 000 xy 2 xx xyxxyxy PHƢƠNG PHÁP SỐ-Bài 6 9 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 Sai số địa phƣơng của thuật toán Taylor bậc k là ...,,),,( ! ),(' ! ),(),( )( 21kyxf k h yxf 2 h yxfyxT 1k 1k k     hxxy 1k h 1k yfh E nn 1k 1kk1k          , )( )!()!( ))(,( )( )( PHƢƠNG PHÁP SỐ-Bài 6 10 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 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 do đó và  thuật toán Taylor bậc cao hơn khó đƣợc chấp nhận 11ykđmãnathy x y x 1 y 2 2  )(/' á ' ' ),(',),( yy2 x y x y x 2 yxfy x y x 1 yxf 23 2 2    '),( f 2 h fyxT2         ),('),( iiiii1i yxf 2 h yxfhyy PHƢƠNG PHÁP SỐ-Bài 6 11 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 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] có sai số địa phƣơng  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 ),( iii1i yxhfyy  )('' y 2 h E 2  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 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 CÁC PHƢƠNG PHÁP SỐ GIẢI BT CÔSI (2) 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 y’ = f(x, y), y(x0) = y0 Ƣớc lƣợng hàm f(x, y) tại các điểm đƣợ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) Sử dụng CT Euler, tính độ nghiêng ở xi+1 yi+l = yi + 1/2(k1 + k2) • Sai số địa phƣơng )()()( 42 yyx 2 yyxyxx 3 1i1i hOffff2ffff2f 12 h yxy   PHƢƠNG PHÁP SỐ-Bài 6 15 CÁC PHƢƠNG PHÁP SỐ GIẢI BT CÔSI (3) 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 CÁC PHƢƠNG PHÁP SỐ GIẢI BT CÔSI (4) 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 CÁC PHƢƠNG PHÁP SỐ GIẢI BT CÔSI (5) 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 4i0)0( và )4,0 (i0)0( )1(i 4 (i) 4    PHƢƠNG PHÁP SỐ-Bài 6 18 CÁC PHƢƠNG PHÁP SỐ GIẢI BT CÔSI (6) 19 trong đó  4321i1i kk2k2k 6 1 yy  )kyh,hf(xk k 2 1 yh, 2 1 xhfk k 2 1 yh, 2 1 xhfk )y,hf(xk 3ii4 2ii3 1ii2 ii1                 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 : PHƢƠNG PHÁP SỐ-Bài 6 CÁC PHƢƠNG PHÁP SỐ GIẢI BT CÔSI (7) 20 yi+1 xi yi xi+1 y xxi+1/2 k1 k4ktổ hợp k3 k2 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 PHƢƠNG PHÁP SỐ-Bài 6 CÁC PHƢƠNG PHÁP SỐ GIẢI BT CÔSI (8) 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 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  y2 = y(1.4) ≈ y(1.2) + Δy1 = 2.291957 PHƢƠNG PHÁP SỐ-Bài 6 21  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 CÁC PHƢƠNG PHÁP SỐ GIẢI BT CÔSI (9) 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 CÁC PHƢƠNG PHÁP SỐ GIẢI BT CÔSI (10) 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      )z,y,x(g'z )z,y,x(f'y 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 CÁC PHƢƠNG PHÁP SỐ GIẢI BT CÔSI (11) k1 = hf(xi, yi, zi) k2= hf(xi+h/2, yi+k1/2, zi+l1/2) k3 = hf(xi+h/2, yi+k2/2, zi+l2/2) k4 = hf(xi+1, yi+k3, zi+l3) l1 = hg(xi, yi, zi) l2= hg(xi+h/2, yi+k1/2, zi+l1/2) l3 = hg(xi+h/2, yi+k2/2, zi+l2/2) l4 = hg(xi+1, yi+k3, zi+l3) 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 Ví dụ: Giải BT Côsi của hệ dƣới đây trên [0; 0,4] với h = 0,1 với các điều kiện đầu y(0) = 0 và z(0) = 1      yzxx'z z'y PHƢƠNG PHÁP SỐ-Bài 6 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 PHƢƠNG PHÁP SỐ-Bài 6 25 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 CÁC PHƢƠNG PHÁP SỐ GIẢI BT CÔSI (13) 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 tƣơng đƣơng với các đ/k đầu y(0) = 0 và z(0) = 1. Giải hệ này tìm được y      yzxx'z z'y PHƢƠNG PHÁP SỐ-Bài 6 26 CÁC PHƢƠNG PHÁP SỐ GIẢI BT CÔSI (14) 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