Bài giảng Phương pháp số trong tính toán cơ khí - Bài 7: Phương trình trình vi phân thường bậc I - Trường Đại học Công nghệ Thành phố Hồ Chí Minh

m riêng 2 2 2 d v  u  u  6tv 1 2  2  0 dt 2 y  x Gồm 1 hoặc nhiều đạo Gồm một hay nhiều đạo hàm hàm toàn phần của của các hàm ẩn số theo một biến độc lập các hàm ẩn số theo nhiều biến độc lập Khái niệm phương trình vi phân 5 - Phương trình vi phân (PTVP) là một phương trình mà chứa đạo hàm của một hàm số chưa biết. Nghiệm của PTVP chính là hàm số mà thỏa mãn PTVP đó. - PTVP mà chỉ có một biến số thì được gọi là Phương trình vi phân thường (Ordinary Differential Equation - ODE) x là biến độc lập; y là biến phụ thuộc  Phương trình có chứa x, y, dy/dx dy f x, y , 0 (1)  PTVP thường bậc 1 là phương trình chứa dx hàm x, hàm y và dy/dx dy - Ví dụ 1: ax2  by  0 PTVP thường bậc 1 (dy/dx), tuyến tính (đối với y) dx dy - Ví dụ 2: a  y  x  b y  0 PTVP thường bậc 1 (dy/dx), phi tuyến (đối với y) dx Phương trình vi phân đầy đủ 6 d - Ví dụ 1: y x  ax2  b  y x  0 dx d - Ví dụ 2: y x  a  y x  x  b y x  0 dx y(x) là hàm ẩn x – biến độc lập Bậc của phương trình vi phân 7 d x t x t et dt PTVP bậc 1 d2 x t d x t 5  2x t  cos t PTVP bậc 2 dt2 dt d3 x t d x t  21xt4    PTVP bậc 3 dt3 dt Phương trình vi phân tuyến tính và phi tuyến 8 Một PTVP là tuyến tính nếu hàm ẩn số và đạo hàm của nó xuất hiện với lũy thừa bằng 1. Không có tích của hàm số và/hoặc đạo hàm của nó. d x t t - Ví dụ 1: x t e - PTVP tuyến tính dt d2 x t d x t - Ví dụ 2: 5  2t2 x t  cos t - PTVP phi tuyến dt2 dt 2 d x t 2 - Ví dụ 3: 2t x t cos t - PTVP phi tuyến dt Phương trình vi phân tuyến tính và phi tuyến 9 Một PTVP là phi tuyến nếu hàm ẩn số và đạo hàm của nó xuất hiện với lũy thừa khác 1. Ngoài ra có thể có tích của hàm số và/hoặc đạo hàm của nó. d x t - Ví dụ 4: cos x t et - PTVP phi tuyến dt d2 x t d x t - Ví dụ 5: 52 xt   - PTVP phi tuyến dt2 dt d2 x t d x t - Ví dụ 6: 51 xt   - PTVP phi tuyến dt2 dt Điều kiện phụ để giải PTVP 10 Nghiệm của phương trình vi phân: d2 x t 4x t  0  x t  c sin 2 t  c cos 2 t dt 2 12 Tất cả các hàm số trên là nghiệm của PTVP, chúng khác nhau ở các hằng số c1, c2. Để xác định chính xác c1, c2 cần có thêm các điều kiện phụ: d2 x t  40xt  - PTVP bậc 2  dt 2 Để giải PTVP bậc n  x0  a ; x 0 b 2 điều kiện phụ ta cần n điều kiện b phụ x t sin 2 t  a cos 2 t 2 Phân loại PTVP theo điều kiện phụ 11 Bài toán giá trị ban đầu Bài toán giá trị biên (Initial Value Problem) (Boundary Value Problem) • Các điều kiện thì không ở 1 Tất cả điều kiện là ở 1 điểm của biến độc lập điểm của biến độc lập • Giải bài toán này khó hơn bài toán giá trị ban đầu x2 x  x  e2t x 2x  x  e2t xx(0) 1, (0) 2.5 x(0) 1, x(2) 1.5 Giống Khác nhau nhau Phương trình vi phân bậc I 12 dy  f x, y với điều kiện ban đầu: y x00  y (2) dx dy f x, y   1.2 y  7 e0.3x dx - Trường hướng - Các họ đường cong Phương trình vi phân bậc I 13 dy 0.3x  f x, y   1.2 y  7 e Tại mỗi điểm có toạ độ (x,y), ta tính được dx f(x,y), và nó chính là độ dốc (Slope) dy/dx 0xy  4; 2   6 của hàm y(x). Độ nghiêng của các véctơ vẽ tại mỗi điểm chính là độ dốc. Tập hợp các véctơ đó trong khoảng x Є [0; 4]; y Є [2; 6] Với mỗi một điều kiện ban đầu ta có một họ đường cong lời giải: y 26  y 03  y 1.5  2 Phương trình vi phân bậc I 14 dy   f x, y dx (2)  y x00  y . Lời giải số của PTVP bậc I (2) là một tập hợp các điểm rời rạc mà xấp xỉ được hàm số y(x). Khi PTVP được giải bằng phương pháp số thì đề bài phải giới hạn miền lời giải. Ví dụ x Є [a; b]. Tùy thuộc vào phương pháp số được dùng để giải phương trình, số điểm giữa khoảng a và b có thể được xác định. . Ví dụ miền lời giải [a; b] có thể được chia ra làm n khoảng nhỏ đều nhau của biến độc lập từ x0 = a cho đến xn = b. . Các lời giải bao gồm các giá trị của biến phụ thuộc được xác định tại mỗi giá trị của biến độc lập .  Lời giải sẽ là tập hợp các điểm (x0, y0), (x1, y1),,(xn, yn) Phân loại phương pháp giải PTVP thường bậc I 15 Phương pháp số cho PTVP Phương pháp đơn bước Phương pháp đa bước (Single-step Method) (Multiple-step Method) Ước tính nghiệm ở một bước Ước tính nghiệm ở một bước cụ thể dựa trên các thông tin cụ thể dựa trên các thông tin của nhiều hơn một bước của một bước trước trước Lời giải số là một quy trình hoặc tính toán ước tính nghiệm của lời giải chính xác tại một tập hợp các điểm rời rạc. Quá trình giải tăng dần theo từng bước. Bắt đầu từ điểm mà giá trị ban đầu được đưa ra. Sau đó sử dụng lời giải tại điểm đầu tiên để tính được lời giải tại điểm thứ hai gần đó. Tiếp theo là lời giải tại điểm thứ ba, thứ tư, v.vÝ tưởng sử dụng giá trị hàm số tại một vài điểm trước có thể đem đến một ước tính tốt hơn cho việc tìm lời giải. Phân loại phương pháp giải PTVP thường bậc I 16 Phương pháp số cho PTVP Phương pháp tường minh Phương pháp ẩn tàng (Explicit Method) (Implicit Method) yi11 F x i,, x i y i  (3) yi1 F x i,,, y i x i  1 y i  1  (4) yi+1 chưa biết và xuất hiện ở cả 2 vế. Nếu F là Trong công thức tường minh, vế hàm tuyến tính của yi+1 thì có thể đưa (4) về (3). phải của phương trình chỉ chứa Còn nếu là phi tuyến thì phải sử dụng kiến thức những giá trị đã biết: giải phương trình phi tuyến để tìm yi+1. Các phương pháp ẩn ngầm đưa đến kết quả chính  xi, xi+1, yi đều đã biết xác hơn so với các phương pháp tường minh nhưng đòi hỏi nhiều nỗ lực tính toán hơn. Phương pháp đơn bước – tường minh 17 dy y  fxyya ,;;    yxaxbx     dx 00 n x1 x 0  h; x 2  x 0  2 h ; ; xn  x 0  nh Trong phương pháp đơn bước – tường minh, Lời giải số ước tính nghiệm của lời giải số (xi+1,yi+1) được tính từ lời giải đã biết (xi,yi). xii1  x h Độ dốc = ϕi  (5) y y   h Lời giải  i1 i i chính xác Trong đó: Đường cong h – cỡ bước của lời giải ϕi=độ dốc – hằng số ước tính giá trị của dy/dx chính xác trong khoảng từ xi đến xi+1. Các phương pháp khác nhau ở chỗ cách tính độ dốc ϕ Phương pháp Euler tường minh 18 dy Độ dốc =  f x, y  (6) idx i i xx i 1.B1: Xác định hàm f(x,y) 2.B2: Nhập số liệu điều kiện ban đầu x0, y0 3.B3: Xác định cỡ bước h, số bước n 4.B4: Cho i từ 0 đến n-1, thực hiện bước 5 và Lời giải 6 chính xác 5.B5: Tính: i f x i, y i  Lời giải số yi1  y i  h  i (7) xii1  x h 6.B6: Xuất kết quả tập hợp (x0,y0), (x1,y1),,(xn,yn)  Chính là lời giải số của PTVP Phương pháp Euler tường minh 19 dy 0.3x  f x, y   1.2 y  7 e dx  yx0  3;0   4. 0.3x g x 1; f10 x   1.2; f x  7 e  fx1   1.2 6 F x  dx  dx  1.2 x   x  gx  15  9 Fx fx  70 x  e  0 dx e1.2x 77 e 0.3 x dx  e 0.9 x dx  e10      gx  9 69   xx 570 10 69 y x  e  C  e 43 xx 70 43  9 C    y x  e5  e 10   9 9 9 y 03  Phương pháp Euler tường minh 20 dy 0.3x 3. i=2:  f x, y   1.2 y  7 e  f x, y   1.2  4.892477918  7 e0.31.0   0.685245956 dx 2 2 2 y y  h  4.892477918  0.5   0.685245956 yx0  3; 0   4 3 2 2   4.54985494 x32 x  h 1.0  0.5  1.5 ba40 4. i=3: hn0.5     8 0.31.5  f x, y   1.2  4.54985494  7 e   0.996428866 h 0.5 3 3 3  y4 y 3  h  3 4.54985494  0.5   0.996428866  4.051640507 1. i=0: 0.3 0 0f x 0, y 0    1.2  3  7 e  3.4 x43 x  h 1.5  0.5  2.0 y y  h  3  0.5  3.4  4.7 1 0 0 5. i=4: x x  h 0  0.5  0.5 10  f x, y   1.2  4.051640507  7 e0.3 2.0   1.020287156 2. i=1: 4 4 4  0.3 0.5 y5 y 4  h  4 4.051640507  0.5   1.020287156 1f x 1, y 1    1.2  4.7  7 e  0.384955835  3.541496929 y2 y 1  h  1 4.7  0.5  0.384955835  4.892477918 x54 x  h 2.0  0.5  2.5 x21 x  h 0.5  0.5  1.0 Phương pháp Euler tường minh 21 0.3 2.5 6. i=5: 5f x 5, y 5    1.2  3.541496929  7 e   0.943230445 y6 y 5  h  5 3.541496929  0.5   0.943230445  3.069881706 x65 x  h 2.5  0.5  3.0  f x, y   1.2  3.069881706  7 e0.3 3.0   0.837870429 7. i=6: 6 6 6  y7 y 6  h  6 3.069881706  0.5   0.837870429  2.650946492 x76 x  h 3.0  0.5  3.5 0.3 3.5 8. i=7: 7f x 7, y 7    1.2  2.650946492  7 e   0.731571546 y8 y 7  h  7 2.650946492  0.5   0.731571546  2.285160719 x87 x  h 3.5  0.5  4.0 Phương pháp Euler tường minh 22 69  xx70 43 y x  e5  e 10  99 Phương pháp Euler tường minh (h=0.5) Lời giải chính xác Phương pháp Euler ẩn tàng 23 dy Độ dốc =  f x, y  i1dx i  1 i  1 (8) xx i1 1.B1: Xác định hàm f(x,y) 2.B2: Nhập số liệu điều kiện ban đầu x0, y0 3.B3: Xác định cỡ bước h, số bước n 4.B4: Cho i từ 0 đến n-1, thực hiện bước 5÷7 Lời giải số 5.B5: Tính: xii1  x h y y  h  f x, y  (9) Lời giải i1 i i  1 i  1 chính xác 6. B6: Giải (9) để tìm yi+1 bằng các pp số 7. B7: Xuất kết quả tập hợp (x0,y0), (x1,y1),,(xn,yn)  Chính là lời giải số của PTVP Phương pháp Euler ẩn tàng 24 dy 0.3x  f x, y   1.2 y  7 e 2. i=1: dx   yx0  3; 0   4 x x  h 0.5  0.5  1.0  21  0.31.0 y2 y 1  h  f x 2, y 2  3.757798698  0.5   1.2  y 2  7 e   3.757798698 3.5e0.31.0 ba40 y   3.969164044 hn0.5     8  2 1.6 h 0.5 1. i=0: 3. i=2:   x x  h 0  0.5  0.5 x x  h 1.0  0.5  1.5  10  32  0.3 0.5  0.31.5 y1 y 0  h  f x 1, y 1  3  0.5   1.2  y 1  7 e  y3 y 2  h  f x 3, y 3  3.969164044  0.5   1.2  y 3  7 e   0.3 0.5  0.31.5  3 3.5e  3.969164044 3.5e y   3.757798698 y3   3.875539109  1 1.6  1.6 6. i=5: 25 Phương pháp Euler ẩn tàng  x x  h 2.5  0.5  3.0  65  0.3 3.0 y6 y 5  h  f x 6, y 6  3.297512707  0.5   1.2  y 6  7 e   0.3 3.0  3.297512707 3.5e 4. i=3: y6   2.950316572  1.6   7. i=6: x43 x  h 1.5  0.5  2.0  y y  h  f x, y  3.875539109  0.5   1.2  y  7 e0.3 2.0   4 3 4 4 4    x76 x  h 3.0  0.5  3.5 3.875539109 3.5e0.3 2.0    0.3 3.5 y4   3.622737397 y y  h  f x, y  2.950316572  0.5   1.2  y  7 e  1.6  7 6 7 7  7   2.950316572 3.5e0.3 3.5 y   2.609436684  7 5. i=4:  1.6 8. i=7:    x54 x  h 2.0  0.5  2.5   x87 x  h 3.5  0.5  4.0  0.3 2.5  y5 y 4  h  f x 5, y 5  3.622737397  0.5   1.2  y 5  7 e   0.3 4.0 y8 y 7  h  f x 8, y 8  2.609436684  0.5   1.2  y 8  7 e   0.3 2.5  3.622737397 3.5e  0.3 4.0 y5   3.297512707  2.609436684 3.5e  1.6 y8   2.289760266  1.6 Phương pháp Euler ẩn tàng 26 69  xx70 43 y x  e5  e 10  99 Phương pháp Euler tường minh (h=0.5) Lời giải chính xác Phương pháp Euler ẩn tàng (h=0.5) Phương pháp Taylor bậc cao 27 Triển khai y(xi+1) trong dãy Taylor quanh điểm xi: 121 1kkkk 1   1 y xi1   y x i  hy x i  h y  x i   h y x i  h y  i  (10) 2!kk !  1 ! Trong đó ξi Є [xi, xi+1], h= xi+1– xi  1 1 1 y x  y x  h y x  hy  x   hkk1 ykk  x  h y 1   i1 i i2! ikk ! i  1 ! i y f x, y yi1 i  ii pk x i, y i  yi1  y i  h  p k x i, y i  (11) Trong đó: 1 1 1 p x,,,,, y  f x y  hf x y  h21 f  x y   hk  fk 1  x y  (12) kiiii2! ii 3! iik ! ii Phương pháp Taylor bậc hai 28 xii1  x h y y  h  p x, y  i1 i k i i (13) 1 1 1 p x,,,,, y  f x y  hf x y  h21 f  x y   hk  fk 1  x y  kiiii2! ii 3! iik ! ii - Phương pháp Euler chính là khi k=1: p1  xi,, y i  f x i y i  i - Phương pháp Taylor bậc hai k=2: 1 pxy2  i,,,,, i  fxy i i  hgxyfxy   i i i i  1  2 p2  x,,, y  f x y hf x y  2  yi12  y i  h  p x i, y i  fxy ,,,,,,  gxyy  gxyfxy     xii1  x h i từ 0 đến n-1 (14) Phương pháp Taylor bậc hai 29 dy 0.3x  0.3x p x, y  0.84 y  4.375 e i  f x, y   1.2 y  7 e 2  i i i dx  y y  h  p x, y   i12 i i i yx0  3; 0   4 xii1  x h hn0.5   8 1. i=0: 0.3 0 p2 x 0, y 0    0.84  3  4.375 e  1.855 0.3x f x, y   1.2 y  2.1 e y1 3  0.5  1.855  3.9275 0.3xx 0.3 x 0  0.5  0.5  1.2  1.2y  7 e  2.1 e 1 0.3x 1.44ye 10.5 2. i=1: 1 0.3 0.5 p2  x,,, y  f x y hf x y p x, y  0.84  3.9275  4.375 e  0.466497396 2 2 1 1  0.5 y2 3.9275  0.5  0.466497396  4.160748698  1.2y  7 e0.3xx  1.44 y  10.5 e 0.3  2 x2 0.5  0.5  1.0  0.84ye  4.375 0.3x Phương pháp Taylor bậc hai 30 0.31.0 3. i=2: p2 x 2, y 2    0.84  4.160748698  4.375 e   0.25394919 y3 4.160748698  0.5   0.25394919  4.033774103 x3 1.0  0.5  1.5 4. i=3: 0.31.5 p2 x 3, y 3    0.84  4.033774103  4.375 e   0.598747083 y4 4.033774103  0.5   0.598747083  3.734400561 x4 1.5  0.5  2.0 5. i=4: 0.3 2.0 p2 x 4, y 4    0.84  3.734400561  4.375 e   0.735845563 y5 3.734400561  0.5   0.735845563  3.366477780 x5 2.0  0.5  2.5 Phương pháp Taylor bậc hai 31 6. i=5: 0.3 2.5 p2 x 5, y 5    0.84  3.366477780  4.375 e   0.761237667 y6 3.366477780  0.5   0.761237667  2.985858946 x6 2.5 0.5 3.0 7. i=6: 0.3 3.0 p2 x 6, y 6    0.84  2.985858946  4.375 e   0.729379253 y7 2.985858946  0.5   0.729379253  2.621169319 x7 3.0  0.5  3.5 8. i=7: 0.3 3.5 p2 x 7, y 7    0.84  2.621169319  4.375 e   0.670804575 y8 2.621169319  0.5   0.670804575  2.285767031 x8 3.5  0.5  4.0 Phương pháp Taylor bậc hai 32 69  xx70 43 y x  e5  e 10  99 Phương pháp Euler tường minh (h=0.5) Lời giải chính xác Phương pháp Taylor bậc 2 (h=0.5) Phương pháp Euler ẩn tàng (h=0.5) Phương pháp Runge-Kutta bậc 2 (RK2) 33 dy y  fxyya ,;;    yxaxbx     dx 00 n x1 x 0  h; x 2  x 0  2 h ; ; xn  x 0  nh Đây là một hệ phương trình 4 ẩn (a , a , b , c ) nhưng chỉ có 3 yi1  y i  h   x i, y i  1 2 1 11 phương trình. Nếu một trong số các giá trị là hằng số thì 3 giá trị (15) yi1 y i  h  a 1 k 1, i  a 2 k 2, i  (RK2) còn lại sẽ có thể được tính. Vì vậy chúng ta có một vài phiên bản của phương pháp RK2: k1,i f x i, y i  (16)  - Phương pháp Euler Cải tiến k f x  b h, y  c k h  2,i i 1 i 11 1, i  - Phương pháp Heun - Phương pháp Ralston 11 (17) a a 1; a b  ; a c  1 2 2 122 2 11 Phương pháp Runge-Kutta bậc 2 (RK2) 34  k f x, y a1  0 1,i i i    17 1 1 1 a1,   b   k  f x  h y  k h 2 1 2,i i i 1, i (18) 2 2 2 1 y y  h  k c   i1 i 2, i  11 2 xii1  x h 11 (19) yi1  yhfx i   i  hy,, i  hfxy   i i  22 Phương pháp Runge-Kutta bậc 2 (RK2) 35 dy f x, y   1.2 y  7 e0.3x xx0.5 dx ii1  0.5 0.5 0.3xi yx0  3; 0   4; y y 0.5   1.2 y    1.2 y  7 e0.3xi  7 e 2 i1 i i i   2 hn0.5   8  x x h ii1 xx0.5 11  ii1 y yhfx    hy,,  hfxy  i1 i i i i i   0.3x 0.25 22 0.3xi  i  yii1 0.58 y  1.05 e  3.5 e 1. i=0: 2. i=1: x1 0  0.5  0.5 x2 0.5  0.5  1.0  0.3 0 0.3 0 0.25  0.3 0.5 0.3 0.5 0.25 y1 0.58  3  1.05 e  3.5 e  3.937102202 y2 0.58  3.937102202  1.05 e  3.5 e    4.174582668 Phương pháp Runge-Kutta bậc 2 (RK2) 36 69  xx70 43 y x  e5  e 10  99 Phương pháp Euler tường minh (h=0.5) Lời giải chính xác Phương pháp Euler cải tiến (h=0.5) Phương pháp Taylor bậc 2 (h=0.5) Phương pháp Euler ẩn tàng (h=0.5) Phương pháp Runge-Kutta bậc 2 (RK2) 37  1  a   1 k1,i f x i, y i  2  1 17   a2  b 1 1,  k 2,i  f x i  h y i  k 1, i h 2 (20) c 1 11 1 yi1 y i  h  k 1, i  k 2, i    2 xii1  x h h (21) yi1  y i   fxy i,,, i  fxhyhfxy i  i    i i  2  Phương pháp Runge-Kutta bậc 2 (RK2) 38  Euler h f xi,, y i  f x i11 y i  y y   fxy ,,,  fxhyhfxy      y  yh  (22) i11 i22 i i i i i i i i xi1 Euler yi1 Phương pháp Runge-Kutta bậc 2 (RK2) 39 xx0.5  ii1  0.5 0.3xxii 0.3 0.3xi 0.5 yi1  y i  1.2 y i  7 e   1.2 y i  0.5  1.2 y i  7 e  7 e  2      xii1  x h h yi1  y i   fxy i,,, i  fxhyhfxy i  i    i i  2  xxii1 0.5  0.3xi 0.3xi 0.5 yii1 0.58 y  0.7 e  1.75 e 1. i=0: 2. i=1: x2 0.5  0.5  1.0 x1 0  0.5  0.5  0.3 0.5 0.3 0.5 0.5  0.3 0 0.3 0 0.5 y0.58  3.946238959  0.7 e  1.75 e y0.58  3  0.7 e  1.75 e  3.946238959  2  1    4.187746066 Phương pháp Runge-Kutta bậc 2 (RK2) 40 69  xx70 43 y x  e5  e 10  99 Phương pháp Euler tường minh (h=0.5) Lời giải chính xác Phương pháp Heun (h=0.5) Phương pháp Euler cải tiến (h=0.5) Phương pháp Taylor bậc 2 (h=0.5) Phương pháp Euler ẩn tàng (h=0.5) Phương pháp Runge-Kutta bậc 2 (RK2) 41  1  a1  k f x, y  3 1,i i i    217  3 3 3 a2  b 1   k 2,i  f x i  h, y i  k 1, i h 3 4 4 4 (22) 3 c  1 11 yi1 y i  h  k 1, i  2 k 2, i   4  3 xii1  x h h 33 (23) yi1  y i   fxy i, i  2 fx i  hy , i  hfxy   i , i  3 4 4 Phương pháp Runge-Kutta bậc 2 (RK2) 42 xxii1 0.5   3 0.3x   0.5  0.5 3 i 0.3xxii 0.3 4 yi1  y i 1.2 y i 7 e  2 1.2 y i  0.5 1.2 y i 7 e  7 e 34   xx0.5  ii1  770.3x 0.3x 0.375 y0.58 y  ei  e i  ii1 60 3 1. i=0: 2. i=1: x2 0.5  0.5  1.0 x1 0  0.5  0.5    770.3 0.5 0.3 0.5 0.375  77 y2 0.58  3.941727143  e  e y0.58  3  e0.3 0  e0.3 0 0.375  3.941727143 60 3  1 60 3    4.181245857 Phương pháp Runge-Kutta bậc 2 (RK2) 43 69  xx70 43 y x  e5  e 10  99 Phương pháp Euler tường minh (h=0.5) Lời giải chính xác Phương pháp Ralston (h=0.5) Phương pháp Heun (h=0.5) Phương pháp Euler cải tiến (h=0.5) Phương pháp Taylor bậc 2 (h=0.5) Phương pháp Euler ẩn tàng (h=0.5) Phương pháp Runge-Kutta bậc 3 (RK3) 44 dy y  fxyya ,;;    yxaxbx     dx 00 n x1 x 0  h; x 2  x 0  2 h ; ; xn  x 0  nh yi1  y i  h   x i, y i  (24) yi1 y i  h  a 1 k 1, i  a 2 k 2, i  a 3 k 3, i  (RK3)  k1,i f x i, y i   (25) k2,i f x i  b 1 h, y i  c 11 k 1, i h  k fx  bhy,  ckhckh   3,i i 2 i 21 1, i 22 2, i  Phương pháp Runge-Kutta bậc 3 (RK3) 45 k1,i f x i, y i    11 k2,i f x i  h, y i  hk 1, i  22  (26) k3,i f x i  h,2 y i  hk 1, i  hk 2, i    h yi1 y i   k 1, i 4 k 2, i  k 3, i   6 xii1  x h h 1 1   1 1  yii1  y f x ii, y  4 f x ii h , y hf x ii , y   f x iiii h , y hf x , y  2 hf  x ii h , y hf x ii , y   6 2 2   2 2  (27) Phương pháp Runge-Kutta bậc 3 (RK3) 46 dy 0.3x  f x, y   1.2 y  7 e dx  y0  3; 0  x  4; h  0.5  n  8 1. i=0: 0.3 0 k1,0 f x 0, y 0    1.2  3  7 e  3.4   1 1   1  0.3 0.25 k2,0 f x 0 h, y 0 hk 1,0   f  0.25,3  0.5 3.4  f 0.25,3.85  1.2 3.85 7 e  1.874204404  2 2   2   k3,0 f x 0 h, y 0 hk 1,0 2 hk 2,0   f 0.5,3  0.5 3.4 2 0.5 1.874204404  f  0.5,3.174204404  0.3 0.5   1.2  3.174204404  7e  2.21591055  h 0.5 y1 y 0 k 14 k 2,0 k 3,0  3 3.4 4 1.874204404  2.21591055  4.092727347  66 Tương tự cho các vòng lặp sau Phương pháp Runge-Kutta bậc 3 (RK3) 47 69  xx70 43 y x  e5  e 10  99 Phương pháp Euler tường minh (h=0.5) Lời giải chính xác Phương pháp RK3 cổ điển (h=0.5) Phương pháp Ralston (h=0.5) Phương pháp Heun (h=0.5) Phương pháp Euler cải tiến (h=0.5) Phương pháp Taylor bậc 2 (h=0.5) Phương pháp Euler ẩn tàng (h=0.5) Phương pháp Runge-Kutta bậc 3 (RK3) 48 k1,i f x i, y i    11 k2,i f x i  h, y i  hk 1, i  33  (28)  22 k3,i f x i  h, y i  hk 2, i  33  h yi1 y i   k 1, i  3 k 3, i   4 xii1  x h h 2 2 1 1 (29) yi1  y i   f x i, y i  3 f x i  h , y i  hf x i  h , y i  hf x i , y i  4 3 3 3 3 Phương pháp Runge-Kutta bậc 3 (RK3) 49 dy 0.3x  f x, y   1.2 y  7 e dx  y0  3; 0  x  4; h  0.5  n  8 1. i=0: k f x, y   1.2  3  7 e0.3 0  3.4  1,0 0 0 1  0.3 1 1   0.5 0.5   1 3.4  3.4 6 k2,0 f x 0 h, y 0 hk 1,0   f  ,3  3.4  f  ,   1.2 7 e  5.978605972  3 3   3 3   6 6  6  2 2   1 1   1 k3,0 f x 0  h, y 0  hk 2,0   f  ,3   5.978605972   f  ,1.992868657  3 3   3 3   3 1  0.3   1.2  1.992868657  7e 3  3.942419538   h 0.5 y1 y 0   k 1,0 3 k 3,0   3   3.4  3  3.942419538  4.093407327  44 Tương tự cho các vòng lặp sau Phương pháp Runge-Kutta bậc 3 (RK3) 50 69  xx70 43 y x  e5  e 10  Phương pháp 99  RK3-Heun (h=0.5) Lời giải chính xác Phương pháp Runge-Kutta bậc 4 (RK4) 51 dy y  fxyya ,;;    yxaxbx     dx 00 n x1 x 0  h; x 2  x 0  2 h ; ; xn  x 0  nh yi1  y i  h   x i, y i  (30) yi1 yhak i   11, i  ak 22, i  ak 33, i  ak 44, i  (RK4) k1,i f x i, y i   k f x  b h, y  c k h  2,i i 1 i 11 1, i  (31)  k3,i fx i  bhy 2, i  ckhckh 21 1, i  22 2, i   k fx  bhy,  ckhckhckh    4,i i 3 i 31 1, i 32 2, i 33 3, i  Phương pháp Runge-Kutta bậc 4 (RK4) 52 k1,i f x i, y i    11 k2,i f x i  h, y i  hk 1, i  22   11 k3,i f x i  h, y i  hk 2, i (32)  22  k4,i f x i  h, y i  hk 3, i   h yi1 y i   k 1, i 22 k 2, i  k 3, i  k 4, i   6  xii1  x h Phương pháp Runge-Kutta bậc 4 (RK4) 53 dy 0.3x  f x, y   1.2 y  7 e dx  y0  3; 0  x  4; h  0.5  n  8 1. i=0: 0.3 0 k1,0 f x 0, y 0    1.2  3  7 e  3.4   11 0.3 0.25 k2,0 f x 0 h, y 0 hk 1,0  f 0.25,3  0.25 3.4 f 0.25,3.85  1.2 3.85 7 e  1.874204404  22   11 k3,0 f x 0  h, y 0  hk 2,0  f 0.25,3  0.25  1.874204404  f  0.25,3.468551101  22    1.2 3.468551101 7e0.3 0.25 2.331943083  k4,0 f x 0  h, y 0  hk 3,0   f 0.5,3  0.5  2.331943083  f  0.5,4.165971542  0.3 0.5   1.2  4.165971542  7e  1.025789985  h 0.5 y1 y 0 k 1,02 k 2,0 2 k 3,0 k 4,0  3 3.4 2 1.874204404  2 2.331943083 1.025789985  66   4.069840413 Tương tự cho các vòng lặp sau Phương pháp Runge-Kutta bậc 3 (RK3) 54 69  xx70 43 y x  e5  e 10  Phương pháp 99  RK4 Cổ điển (h=0.5) Lời giải chính xác Sử dụng lệnh có sẵn của MATLAB 55 69  xx Lời giải 570 10 43 dy 0.3x y x  e  e   f x, y   1.2 y  7 e chính xác 99 dx  y0  3; 0  x  4; h  0.5  n  8 1) Tạo file hàm số xác định hàm dy/dx = f(x,y): function dydx = DiffEq(x,y) dydx = -1.2*y + 7*exp(-0.3*x); end - Giải bằng một trong số các câu lệnh có sẵn của 2) Tạo file script lời giải: MATLAB, ví dụ như ode45 - Giải bằng lời giải chính xác - Nhập vào MATLAB lời giải bằng tay mà ta đã tính dựa trên 1 trong số 9 phương pháp đã học - So sánh kết quả giữa lời giải chính xác và MATLAB - So sánh kết quả giữa lời giải chính xác và lời giải tay 04x format long;  56 h  0.5 tspan = [0:0.5:4]; % Mien cua bien doc lap:Tu a den b voi buoc h y 03  y0 = 3; % Gia tri ban dau cua ham phu thuoc y 69 [x,yMATLAB] = ode45(@DiffEq,tspan,y0) % Su dung ham ode45  xx70 43 y x  e5  e 10  yExact = exp((-6/5)*x).*((70/9)*exp((9/10)*x) - 43/9) % Loi giai 99 chinh xac % Nhap ket qua tinh tay theo mot trong so 9 phuong phap, vi du PP Euler tuong minh y_ExEuler = [3;4.7;4.892477918;4.549854940;4.051640506;3.541496929;3.06988170 6;2.650946492;2.285160719] % Tinh sai so giua loi giai chinh xac va MATLAB error_Exact_MATLAB = abs((yExact - yMATLAB)./yExact)*100 % Tinh sai so giua loi giai chinh xac va mot trong so 9 pp ma ta giai tay, vi du PP Euler tuong minh error_Exact_ExEuler = abs((yExact - y_ExEuler)./yExact)*100 % Ve do thi ham loi giai chinh xac x = [0:0.01:4]; %Chia nho lai khoang x de ve do thi cho muot yExact = exp((-6/5)*x).*((70/9)*exp((9/10)*x) - 43/9); plot(x,yExact,'b') hold on scatter(tspan,yMATLAB,'k') %ve cac diem loi giai cua MATLAB hold on scatter(tspan,y_ExEuler,'d') %Ve cac diem loi giai cua tinh tay 57 Véctơ x Véctơ lời giải y theo câu Véctơ lời giải chính xác lệnh ode45 của MATLAB Phương pháp Euler tường minh (h=0.5) Véctơ lời giải theo 1 Sai số giữa lời giải Sai số giữa lời giải chính xác và lời giải chính xác và lời giải trong số 9 phương Lời giải chính xác pháp số tính bằng tay bằng câu lệnh của tính bằng tay theo 1 MATLAB trong số 9 pp Lời giải theo câu lệnh ode45 của MATLAB Bài toán kỹ thuật 58 Cho ống nước chảy vào bể nước hình trụ có bán kính Rtank. Cùng ngay tại thời điểm đó thì ở dưới đáy bể, nước cũng được thoát ra bởi một vòi ống tròn có bán kính rhole. Cho biết quy luật dòng chảy nước vào bể có dạng: dm in K Kcos Ct (33) dt 12 Cho biết ở thời điểm ban đầu t0=0 thì nước trong bể đang ở độ cao h0. Hãy tính độ cao của bể nước ở thời điểm t* (s) kể từ khi bắt đầu. Cho biết khối lượng riêng của nước là ρ. Để giản lược tính toán, ta chấp nhận giản lược sau của độ cao h khi tính dòng chảy ra khỏi bể: hh 0.008 (34) Giải phương trình vi phân bằng phương pháp giải tích chính xác và một trong số 9 phương pháp số có trong bài học để so sánh. Cho bước lặp bằng 0.5. Khảo sát thời gian t Є [0;10]. 59 - Dòng chảy ra ở vòi ống dưới đáy bể: 2 Arhole  hole s dmout ds - Khối lượng của nước đi qua ống độ dài s: m  A  s    A  out holedt hole dt (34) v - Vận tốc của nước dưới đáy bể khi nước có độ cao h: v 2 gh (35) dm 22 - Dòng chảy ra ở vòi dưới đáy bể là: out   r 2 gh     r  2 g  0.008 h dt hole hole (36) 60 - Độ biến thiên khối lượng nước trong bể: - Khối lượng của nước trong bể khi có độ cao h: m  Atank  h (37) - Độ biến thiên nước trong bể: dm dh  R2  (38) dttank dt - Theo định luật bảo toàn khối lượng: dmdm dm dh in out  R22  K Kcos C  t   r 2 g 0.008 h dt dt dttank dt 12 hole 2 dh rhole KK21  f h, t  2  0.008 2 g  h  2 cos C  t  2   dt Rtank  R tank   R tank (39)  h t00  h dh 100cos 3t 100 61  f h, t   0.3543557534E-5  h     dt 997 997  h00  rhole  0.01 m ; Rtank 1 m ; K1 100 kg s ; 1) Giải bằng tay bằng phương pháp giải tích. Tìm hàm lời giải chính xác (Slide 19) K2 100 kg s ; 2) Tùy từng đề yêu cầu, giải một trong số 9 phương pháp (Bấm 3 máy tính CASIO chứ không được dùng MATLAB hay Excel hỗ   997 kg m ; trợ vì để luyện tập khi đi thi cuối kz). Chú {: bước lặp bằng 1 0.5. Khảo sát thời gian t Є [0;10]  20 điểm C  3 s ; 3) Lập bảng tính tay so sánh lời giải chính xác và lời giải số theo t  0 (s); mẫu (Slide 22/26/32/36/40/43/47/50/54): 0 h0  0 (m);  t  6 (s); 2 4) Sử dụng Code MATLAB như theo mẫu ở Slide 56 để kiểm g  9.81 m s . nghiệm lại và in kết quả, vẽ đồ thị so sánh kết quả (Slide 57). Phương trình vi phân thường bậc I 62 Phương pháp Euler tường minh (h=0.5) Lời giải chính xác Phương pháp RK3 cổ điển (h=0.5) Phương pháp Ralston (h=0.5) Phương pháp Heun (h=0.5) Phương pháp Euler cải tiến (h=0.5) Phương pháp Taylor bậc 2 (h=0.5) Phương pháp Euler ẩn tàng (h=0.5)