Các phương pháp số cho phương trình vi phân

Chia sẻ: Nguyễn Đăng Khoa | Ngày: | Loại File: PDF | Số trang:43

Thêm vào BST

Báo xấu

356
lượt xem 65
download

Download Vui lòng tải xuống để xem tài liệu đầy đủ

Tài liệu tham khảo về Các phương pháp số cho phương trình vi phân...

Chủ đề:

Bình luận(0) Đăng nhập để gửi bình luận!

Lưu

Nội dung Text: Các phương pháp số cho phương trình vi phân

1 M cl c M cl c 1 0.1 Lý do ch n đ tài . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 0.2 M c đích yêu c u . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 0.3 Phương pháp nghiên c u . . . . . . . . . . . . . . . . . . . . . . . . . 3 0.4 Gi i h n c a đ tài . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 1 Phương pháp s cho phương trình vi phân c p 1 4 1.1 Phương pháp Euler đ gi i phương trình vi phân c p 1 . . . . . . . . 4 1.2 Phương pháp Euler c i ti n (Phương pháp Heun) . . . . . . . . . . . 7 1.3 Phương pháp Runge - Kutta . . . . . . . . . . . . . . . . . . . . . . . 9 1.4 Sai s và s đi u ch nh bư c nhãy RKF (Rung - Kutta - Fehlberg) . . 11 2 Phương pháp nhi u nút 14 2.1 Phương pháp Adams - Bashforth . . . . . . . . . . . . . . . . . . . . 14 2.2 Các phương pháp Adams - Moulton . . . . . . . . . . . . . . . . . . . 15 3 Phương pháp cho h phương trình và phương trình vi phân b c cao 17 3.1 Phương pháp Euler cho h phương trình . . . . . . . . . . . . . . . . 17 3.2 Các phương pháp Runge - Kutta cho h phương trình . . . . . . . . . 19 3.3 Phương pháp Runge - Kutta - Nystr¨m (phương pháp RKN) . . . . . 20 o 4 Phương pháp cho các phương trình vi phân đ o hàm riêng 22 4.1 Các phương pháp cho các phương trình vi phân ki u Elip . . . . . . . 23 4.1.1 Các phương trình vi phân cho phương trình Laplace và phương trình Poisson . . . . . . . . . . . . . . . . . . . . . . . . . . . 23 4.1.2 Bài toán Dirichlet . . . . . . . . . . . . . . . . . . . . . . . . . 25 4.1.3 Phương pháp ADI . . . . . . . . . . . . . . . . . . . . . . . . 28 4.2 Bài toán Neumann và nhi u lo n. Biên không chính quy . . . . . . . 30 4.3 Phương pháp cho phương trình Parabol . . . . . . . . . . . . . . . . . 34
2 4.4 Phương pháp cho phương trình Hypebol . . . . . . . . . . . . . . . . 40 5 Ph n k t lu n 42
3 0.1 Lý do ch n đ tài Phương pháp s là m t môn h c b t bu c cho t t c h c viên Cao h c v t lý. Nó cung c p cho h c viên các phương pháp nh m gi i quy t nh ng v n đ liên quan đ n tính toán, các thu t toán và cách gi i g n đúng hay x p x các phương trình, tích phân, vi phân, đ o hàm,... Phương pháp s có th dùng đ gi i g n đúng các bài toán mà các phương pháp thông thư ng không th gi i đư c. H u h t, chúng mã hoá và đưa ra các thu t gi i đ thu n ti n cho vi c l p trình trên máy tính. Đ tìm hi u sâu hơn v môn h c nhóm chúng tôi đư c s phân công c a th y giáo hư ng d n làm đ tài : Các phương pháp s cho phương trình vi phân (numerical methods for differential equations). Đ làm ti u lu n cho môn h c. Do kinh nghi m còn h n ch nên trong ti u lu n không tránh kh i các sai sót kính mong s góp ý chân thành c a quý th y cô và b n đ c. 0.2 M c đích yêu c u Đưa ra các phương pháp cho phương trình vi phân c p m t, cho h phương trình vi phân c p 1, phương trình vi phân c p cao, phương trình vi phân đ o hàm riêng. 0.3 Phương pháp nghiên c u Ch y u là nghiên c u tài li u và th o lu n nhóm đ đưa ra các k t lu n đưa vào ti u lu n. 0.4 Gi i h n c a đ tài Trong ti u lu n này chúng tôi ch trình bày, đưa ra phương pháp và m i phương pháp ch trình bày m t ví d đ minh ho .
4 Chương 1 PHƯƠNG PHÁP S CHO PHƯƠNG TRÌNH VI PHÂN C P 1 1.1 Phương pháp Euler đ gi i phương trình vi phân c p 1 Ta đã bi t phương trình vi phân c p 1 có d ng F (x, y, y ) = 0, và nó có th vi t tư ng minh y = f (x, y ). Thư ng phương trình vi phân ngư i ta cho thêm m t đi u ki n ban đ u v nghi m tho mãn. Trong ph n này ta s xét phương trình vi phân c p 1 và đi u ki n có d ng: y = f (x, y ); y (x0 ) = y0 (1.1) gi s f là bài toán có duy nh t nghi m trên đo n có ch a x0 . Ta s bàn v các phương pháp đ tính các giá tr s c a các nghi m, đây ta ch c n nêu m t cách th c cho nghi m c a các phương trình không kh d ng ho c quá ph c t p đ tính toán. Các phương pháp mà ta th c hi n đ làm đi u này là phương pháp g n đúng theo t ng b ơc (step-by-step motheds). Ta b t đ u t vi c có giá tr y0 = y (x0 ) và ti n hành theo các bư c, tính các giá tr g n đúng c a nghi m y (x) t i các nút "mesh points " x1 = x0 + h ; x2 = x0 + 2 h ; x3 = x0 + 3h; ... v i h là bư c nhãy (step size) c đ nh, ch ng h n 0.1 hay 0.2 tuỳ chúng ta ch n. Tính toán trong m i bư c đư c th c hi n theo cách th c như nhau. Cách th c đó là d a vào vi c khai tri n chu i Taylor: h2 y (x + h) = y (x) + hy (x) + y (x) + · · · (1.2) 2 Xét v i các giá tr nh c a h, các s h ng b c cao h2 , h3 , · · · là r t nh ,ta có th b qua, do v y: y (x + h) ≈ y (x) + hy (x) = y (x) + hf (x, y )
5 (v ph i đ t đư c trên là t phương trình vi phân y (x) = f (x, y )) và theo các quá trình l p. Trong bư c đ u tiên ta tính y1 = y0 + hf (x0 , y0 ) v i cách l y g n đúng y (x1 ) = y (x0 + h). Trong bư c th hai ta tính y2 = y1 + hf (x1 , y1 ) và m t cách g n đúng thì y (x2 ) = y (x0 + 2h) và t ng quát: yn+1 = yn + hf (xn , yn ); (n = 0, 1, 2, · · · ) (1.3) Đây đư c g i là Phương pháp Euler hay Phương pháp Euler-Cauchy. Phương pháp thô sơ này h u như ít đư c dùng trong th c hành nhưng t phương pháp đơn gi n này, ta gi i thích cho nguyên t c dùng các phương pháp d a trên cơ s khai tri n chu i Taylor. Hình 1: Phương pháp Euler 1 Công th c Taylor v i khai tri n d ng y (x + h) = y (x) + hy (x) + 2 h2 y (ξ ) v i x ≤ ξ ≤ x + h. Nó di n t r ng Phương pháp Euler là làm ng n l i các sai s trong m i bư c hay nó làm ng n sai s b ph n (local truncation error) t l v i h2 , đư c vi t O(h2 ). Trên m t kho ng x c đ nh mà ta mu n gi i m t phương trình s các bư c t l v i 1/h. Sai s b ph n hay sai s toàn ph n đây t l v i h2 (1/h) = h1 . Do nguyên nhân này mà phương pháp Euler đư c g i là Phương pháp c p 1. Ngoài ra còn có sai s do làm tròn và do các phương pháp khác, nó gây nh hư ng đ n đ chính xác c a các giá tr y1 , y2 , · · · Ví d 1.1. Dùng phương pháp Euler đ gi i bài toán đi u ki n đ u, ch n h = 0.2 và tính y1 , · · · , y5 y = x + y, y (0) = 0 (1.4)
6 Gi i: f (x, y ) = x + y , và ta th y r ng (1.3) tr thành: yn+1 = yn + 0.2(xn + yn ) B ng 1.1 di n t cách tính toán, giá tr chính xác c a nghi m c a (1.4) là y (x) = ex − x − 1. Trong th c hành tính toán chính xác nghi m là chưa bi t, nhưng s chính xác c a các giá tr có th tìm đư c b ng vi c áp d ng phương pháp Euler t m t l n khác v i bư c nh y 2h = 0.4 và so sánh v i sai s tương ng. K t qu tính này là: xn yn 0.4(xn + yn yn trong b ng (1.1) S sai khác 0.0 0.000 0.000 0.000 0.000 0.4 0.000 0.160 0.040 0.040 0.8 0.160 0.274 0.114 đây sai s là c a h2 , trong mi n t h đ n 2h nó đư c nhân v i 22 = 4, nhưng vì ta ch c n m t n a c a bư c nhãy trư c đó, nên nó ch nhân v i 4/2 = 2. Do đó hi u s 2 2 − 2 = 0.040, sai s 2 c a y2 trong b ng 1.1 (v i k t qu th c là 0.052) và 0.114 c a y4 (chính xác là 0.152). n xn yn 0.2(xn + yn ) Giá tr đúng Sai s 0 0.0 0.000 0.000 0.000 0.000 1 0.2 0.000 0.040 0.021 0.021 2 0.4 0.040 0.088 0.092 0.052 3 0.6 0.128 0.146 0.222 0.094 4 0.8 0.274 0.215 0.426 0.152 5 1.0 0.489 0.718 0.229 B ng 1.1. Phương pháp Euler áp d ng cho (1.4) trong ví d 1.1 và sai s . L a ch n bi n bư c nh y m t cách ng u nhiên trong các mã hi n đ i. Nh ng mã này l a ch n ki u các bi n nhãy hn sao cho sai s c a nghi m không l n hơn giá tr c c đ i TOL (dung h n đư c đ ngh ). V i phương pháp Euler, khi bư c nhãy 1 là hn , sai s b ph n t i xn là c 2 h2 | y (xn ) |. Ta quy đ nh giá tr này b ng dung n h n TOL 12 2T OL h | y (xn ) |= T OL, nên hn = (1.5) 2n | y (xn ) | y (x) ph i khác không trên đo n J : x0 ≤ x ≤ xN . Ch n K là giá tr c c ti u c a | y (xn ) | trên J và gi s K > 0. Giá tr c c ti u | y (xn ) | tương ng v i giá tr √ √ c c đ i h = H = 2TK c a (1.5). Vì v y, 2T OL = H K . Đưa giá tr này vào OL (1.5) ta đ t đư c K hn = ϕ(xn )H, ϕ(xn ) = (1.6) | y (xn ) |
7 Đ i v i các phương pháp khác,l a ch n bi n bư c nh y m t cách ng u nhiên là như nguyên t c này. 1.2 Phương pháp Euler c i ti n (Phương pháp Heun) Khi ta chú ý nhi u đ n các s trong (1.2) ta đ t đư c các phương pháp s b c cao và chính xác. Nhưng đây là v n đ th c hành. N u ta th y = f (x, y (x)) vào (1.2) ta đư c: 1 1 y (x + h) = y (x) + hf + h2 f + h3 f + · · · (1.7) 2 6 đây, vì f ph thu c vào x nên f = fx + fy y = fx + fy f và các đ o hàm c p cao f , f tr nên c ng k nh hơn. Phương pháp bây gi là hu b các tính toán đó mà thay vào đó là tính f b i m t hay nhi u giá tr ph tr đư c ch n thích h p c a (x, y ). "Thích h p" có nghĩa là chúng đư c ch n làm c p c a phương pháp cao có th th c hi n đư c (có đ chính xác cao). Ta nói v hai phương pháp th c hành quan tr ng. Phương pháp th nh t g i là Phương pháp Euler c i ti n hay Phương pháp Euler - Cauchy c i ti n (còn g i là Phương pháp Heun). Trong m i bư c c a phương pháp này ta tính giá tr ph đ u tiên (1.8a) và r i tính giá tr m i (1.8b) ∗ yn+1 = yn + f (xn , yn ) (1.8a) 1 ∗ yn+1 = yn + h[f (xn , yn ) + f (xn+1 , yn+1 )] (1.8b) 2 Hình 2: Phương pháp Euler c i ti n
8 Phương pháp này có cách gi i thích hình h c đơn giãn. Th c t ta có th trong kho ng t xn đ n xn + 1 h ta có nghi m g n đúng y b ng đư ng th ng qua đi m 2 (xn , yn ) v i h s góc f (xn , yn ), và chúng ta ti p t c d c theo đư ng th ng v i h ∗ s góc f (xn+1 , yn+1 ) cho đ n x là xn+1 . Phương pháp Euler - Cauchy c i ti n là m t Phương pháp predictor - corrector, b i vì trong m i bư c ta đ u tiên ta tiên đoán m t giá tr trong (1.8a) và hi u chính cho chính xác nó trong (1.8b). Ta dùng kí hi u k1 = hf (xn , yn ) trong ∗ (1.8a) và k2 = hf (xn+1 , yn+1 ) trong (1.8b) ta có th vi t sơ đ d ng này trong b ng 1.2. Ví d 1.2. V n d ng phương pháp Euler c i ti n đ gi i bài toán (1.4), ch n h = 0.2 như trư c. Gi i: V i bài toán lúc này. k1 = 0.2(xn + yn ) k2 = 0.2(xn + 0.2)yn + 0.2(xn + yn ) 0.2 yn+1 = yn + (2.2xn + 2.2yn + 0.2) + 0.22(xn + yn ) + 0.02. 2 B ng 1.3 cho th y r ng k t qu các k t qu hi n t i là chính xác hơn các k t qu trong ví d 1.1, xem b ng 1.6 . Thu t toán Euler (f, x0 , y0 , h, N ) Thu t toán này tìm nghi m bài toán v i giá tr đ uy = f (x, y ), y (x0 ) = y0 t i các đi m cách đ u x1 = x0 + h, x2 = x0 + 2h, · · · , xN = x0 + N h, đây, f như m t bài toán có duy nh t nghi m trên đo n [x0 , xN ]. Nh p vào: giá tr ban đ u x0 , y0 , bư c đi(bư c nh y) h s bư c đi trong kho ng kh o sát N Giá tr ra: Giá tr g n đúng yn+1 c a nghi m y (xn+1 ) t i xn+1 = x0 +(n +1)h. V i n = 0, 1, · · · , N − 1. For n = 0, 1, · · · , N − 1 do: xn+1 = xn + h k1 = hf (xn , yn ) k2 = hf (xn+1 , yn + k1 ) 1 yn+1 = (k1 + k2 ) 2 L y ra: xn+1 , yn+1 K t thúc T m d ng K t thúc Phương pháp Euler B ng 1.2. Phương pháp Euler c i ti n(Phương pháp Heun)
9 n xn yn 0.22(xn + yn ) + 0.02 Giá tr chính xác Sai s 0 0.0 0.0000 0.0200 0.0000 0.0000 1 0.2 0.0200 0.0684 0.0214 0.0014 2 0.4 0.0884 0.1274 0.0918 0.0034 3 0.6 0.2158 0.1995 0.2221 0.0063 4 0.8 0.4153 0.2874 0.4255 0.0102 5 1.0 0.7027 0.7183 0.0156 B ng 1.3. Phương pháp Euler c i ti n áp d ng cho (1.4) và sai s Sai s b ph n: Sai s b ph n c a phương pháp Euler c i ti n là s h ng h3 . Th t v y, đ t fn = f (xn , y (xn )) và dùng (1.7), ta đư c 1 1 y (xn + h) − y (xn ) = hfn + h2 fn + h3 fn + · · · (1.9) 2 2 L y g n đúng bi u th c trong ngo c c a (1.8b) b ng fn + fn+1 và dùng khai tri n Taylor ta đư c t (1.8b) 1 yn+1 − yn ≈ h[fn + fn+1 ] 2 1 1 = h[fn + (fn + hfn + h2 fn + · · · )] (1.10) 2 2 12 13 = hfn + h fn + h fn + · · · 2 4 Tr (1.9) cho (1.10) ta có sai s h3 h3 h3 fn − fn + · · · = − fn + · · · 6 4 12 h3 1 = h2 . Vì s các bư c nhãy trên đo n x t l v i h , nên sai s toàn ph n là đây h phương pháp Euler c i ti n đư c g i là Phương pháp c p hai 1.3 Phương pháp Runge - Kutta V n còn phương pháp th c hành quan tr ng chính xác hơn n a là Phương pháp b n s c đi n Runge - Kutta, thư ng g i m t cách ng n g n là Phương pháp Hunge - Kutta. Nó đư c di n t trong b ng 1.4. Ta th y r ng trong m i bư c vi c đ u tiên ta tính b n s k1 , k2 , k3 , k4 và ti p đó tính giá tr m i yn+1 . Công th c có v ph c t p d ng đ u tiên, nhưng th c t nó l i d cho m t bài toán. Vi c tính tay f (x, y ) b ng tay thì r t ph c t p nhưng nó s không thành v n đ đ i v i máy tính. Phương pháp này là hoàn h o cho máy tính b i vì không c n đ n các tính toán ban đ u, nó làm sáng t s tích lũy, và dùng l p l i các k t qu m t cách tr c ti p. Nó là s n đ nh. Chú ý, n u f ch ph thu c vào x, ta dùng công th c tích phân Simpson.
10 Ví d 1.3. áp d ng phương pháp Runge - Kutta đ i v i bài toán (1.4) trong ví d 1.1, ch n h = 0.2 như trư c và tính năm bư c. Gi i: V i bài toán này, ta có f (x, y ) = x + y nên ta có k1 = 0.2(xn + yn ), k2 = 0.2(xn + 0.1 + yn + 0.5k1 ) k3 = 0.2(xn + 0.1 + yn + 0.5k2 ), k4 = 0.2(xn + 0.2 + yn + k3 ) T các bi u th c đơn gi n, ta tìm giá tr thích h p c a k1 và thay vào k2 ta đ t đư c k t qu k2 = 0.22(xn + yn ) + 0.02, thay giá tr này vào k3 ta tìm đư c k3 = 0.222(xn + yn ) + 0.022, và cu i cùng thay giá tr này vào k4 ta đư c k4 = 0.2444(xn + yn ) + 0.0444. N u ta dùng các s trên, công th c yn+1 trong b ng 1.4 tr thành yn+1 = yn + 0.2214(xn + yn ) + 0.0214. (1.11) T t nhiên, các quá trình thay vào này là không tiêu bi u cho phương pháp Kunge - Kutta và s không là v n đ chung. B ng 1.5 di n t quá trình tính toán. T b ng 1.6 ta th y r ng các giá tr là chính xác hơn các giá tr trong ví d 1.1 và ví d 1.2 Thu t toán Runge - Kutta (f, x0 , y0 , h, N ) Thu t toán này tìm nghi m bài toán v i giá tr đ uy = f (x, y ), y (x0 ) = y0 t i các đi m cách đ u x1 = x0 + h, x2 = x0 + 2h, · · · , xN = x0 + N h, đây, f như m t bài toán có duy nh t nghi m trên đo n [x0 , xN ]. Nh p vào: giá tr ban đ u x0 , y0 , bư c đi(bư c nh y) h s bư c đi trong kho ng kh o sát N Giá tr ra: Giá tr g n đúng yn+1 c a nghi m y (xn+1 ) t i xn+1 = x0 +(n +1)h. V i n = 0, 1, · · · , N − 1. For n = 0, 1, · · · , N − 1 do: k1 = hf (xn , yn ) 1 1 k2 = hf (xn + h, yn + k1 ) 2 2 1 1 k3 = hf (xn + h, yn + k2 ) 2 2 k4 = hf (xn + h, yn + k3 ) xn+1 = xn + h 1 yn+1 = yn + (k1 + 2k2 + 2k3 + k4 ) 6 L y ra: xn+1 , yn+1 K t thúc T m d ng K t thúc Phương pháp Runge - Kutta B ng 1.4. Phương pháp Runge - Kutta c đi n c p b n
11 106 × sai s n xn yn 0.2214(xn + yn ) + 0.0214 Giá tr đúng y = c ay ex − x − 1 0 0.0 0.000 000 0.021 400 0.000 000 0 1 0.2 0.021 400 0.070 418 0.021 403 3 2 0.4 0.091 818 0.130 289 0.091 825 7 3 0.6 0.222 106 0.203 414 0.222 119 11 4 0.8 0.425 521 0.292 730 0.425 541 20 5 1.0 0.718 251 0.718 282 31 B ng 1.5. Phương pháp Runge - Kutta áp d ng (1.4); tính toán b ng dùng bi u th c (1.11). y = ex − x − 1 x Sai s Phương pháp Euler Phương pháp Phương pháp Euler c i ti n Runge - Kutta 0.2 0.021 403 0.021 0.0014 0.000 003 0.4 0.091 825 0.052 0.0034 0.000 007 0.6 0.222 119 0.094 0.0063 0.000 011 0.8 0.425 541 0.152 0.0102 0.000 020 1.0 0.718 282 0.229 0.0156 0.000 031 B ng 1.6. So sánh s chính xác c a ba phương pháp xét trong bài toán (1.4) v i h = 0.2. 1.4 Sai s và s đi u ch nh bư c nhãy RKF (Runge - Kutta - Fehlberg) ý tư ng này dùng tích phân thích h p cho Rung-Kutta và các phương pháp khác. Trong b ng 1.4 cho RK (Runge - Kutta), n u ta tính m i bư c v i bư c nhãy h và 2h, thì ít lâu sau sai s cho m i bư c b ng 25 = 32 nhân v i giá tr ban đ u. Tuy nhiên, ta ch có m t n a bư c cho 2h, nên s chính xác là 25 /2 = 16. đây sai s x p xĩ v i bư c nhãy h b ng c 1/15 nhân v i đ chênh l ch δ = y − y c a các giá tr tương ng v i c bư c nhãy l n lư t là h và 2h 1 ≈ (y − y ) (1.12) 12 B ng 1.7 minh ho (1.12) v i phương trình vi phân có giá tr đ u y = (y − x − 1)2 + 2, y (0) = 1; (1.13)
12 y y Sai s Chính xác hoá Nghi m chính x (c h) (c 2h) ơc tính (1.12) sai s chính xác 0.0 1.000 000 000 1.000 000 000 0.000 000 000 0.000 000 000 1.000 000 000 0.1 1.200 334 589 0.000 000 083 1.200 334 672 0.2 1.402 709 878 1.402 707 408 0.000 000 765 0.000 000 157 1.402 710 036 0.3 1.609 336 039 0.000 000 210 1.609 336 250 0.4 1.822 792 993 1.822 788 993 0.000 000 267 0.000 000 226 1.822 793 219 B ng 1.7. Phương pháp Runge - Kutta đư c áp d ng vào phương trình vi phân (1.13) và sái s ư c tính (1.12). Nghi m chính xác y = tan x + x + 1 c bư c h = 0.1 và 0 ≤ x ≤ 4. Ta th y r ng sai s ư c tính g n b ng sai s chính xác. Đây là phương pháp cho sai s ư c tính là đơn giãn nhưng không n đ nh. RKF. E. Fhlberg đ xu t và phát tri n sai s đi u chính b ng cách dùng hai phương pháp RK v i các s khác nhau đi t (xn , yn ) đ n (xn+1 , yn+1 ). S chênh l ch c a vi c tính giá tr y t i xn+1 cho sai s ư c tính đư c dùng cho đi u ch nh c bư c nhãy h. E. Fhlberg đưa ra hai công th c RK v i công th c này ch c n 6 giá tr hàm cho m i bư c. Ta trình bày d ng các hàm đây vì RKF đã tr nên quá ph bi n. Ch ng h n, Maple dùng nó (cho h th ng các phương trình vi phân). Phương pháp RK b c năm c a Fhlberg là yn+1 = yn + γ1 k1 + · · · + γ6 k6 (1.14a) 16 6656 28561 9 2 0 − 50 γ= (1.14b) 135 12825 56430 55 γ1 · · · γ6 V i h ng s vectơ γ = Phương pháp RK b c b n c a ông là ∗ ∗ ∗ yn+1 = yn + γ1 k1 + · · · + γ5 k6 (1.15a) ∗ 25 1408 2197 −1 0 γ= (1.15b) 216 2565 4104 5 Trong c hai công th c ta ch dùng 6 giá tr hàm khác nhau k1 = hf (xn , yn ) 1 1 k2 = hf (xn + h, yn + k1 4 4 3 3 9 k3 = hf (xn + h, yn + k1 + k2 8 32 32 12 1932 7200 7296 k4 = hf (xn + h, yn + k1 − k2 + k3 ) (1.16) 13 2197 2197 2197 439 3680 845 k5 = hf (xn + h, yn + k1 − 8k2 + k3 − k4 ) 216 513 4104 1 8 3544 1859 11 k6 = hf (xn + h, yn − k1 + 2k2 − k3 + k4 − k5 2 27 2565 4104 40
13 S sai khác c a (1.14a) và (1.15a) cho sai s ư c tính 1 128 2197 1 2 ∗ ≈ yn+1 − yn+1 = k1 − k3 − k4 + k5 + k6 (1.17) n+1 360 4275 75240 50 55 Ví d 1.4. Runge - Kutta - Fehlberg cho bài toán phương trình vi phân (1.13) ta đ t đư c t (1.14a), (1.14b), (1.15a), (1.15b), (1.16) v i h = 0.1 trong bư c đ u tiên v i 12 s th p phân là k1 = 0.200000000000 k2 = 0.200062500000 k3 = 0.200140756867 k4 = 0.200856926154 k5 = 0.201006676700 k6 = 0.200250418651 ∗ y1 = 1.20033467253 y1 = 1.20033466949 1 = 0.00000000304 Giá tr chính xác đ n 12 s th p phân là y (0.1) = 1.20033467209. Do đó sai s chính xác c a y1 là −4.4 · 10−10 , nh hơn trong b ng 1.7 b ng h ng s 200. B ng 1.8 t ng k t các k t qu cơ b n c a các phương pháp trong chương 1 này. Nó có th các phương pháp này là các sai s n đ nh. Chúng là các phương pháp m t bư c b i vì m i bư c ta dùng d ki n c a bư c trư c m t cách có th t. Phương pháp Hàm giá tr cho m i bươc Sai s toàn ph n Sai s c c b O(h2 ) Euler 1 O(h) O(h2 ) O(h3 ) Euler c i ti n 2 O(h4 ) O(h5 ) RK(b c b n) 4 O(h5 ) O(h6 ) RKF 6 B ng 1.8. Các phương pháp đư c xem xét và b c c a chúng (=sai s toàn ph n c a chúng) 106 × sai s B t đ u Tiên đoán Hi u ch nh Giá tr ∗ n xn yn yn yn chính xác c a yn 0 0.0 0.000 000 0.000 000 0 1 0.2 0.021 400 0.021 403 3 2 0.4 0.091 818 0.091 825 7 3 0.6 0.222 107 0.222 119 12 4 0.8 0.425 361 0.425 529 0.425 541 12 5 1.0 0.718 066 0.718 270 0.718 282 12 6 1.2 1.119 855 1.120 106 1.120 117 11 7 1.4 1.654 885 1.655 191 1.655 200 9 8 1.6 2.352 653 2.353 026 2.353 032 6 9 1.8 3.249 190 2.249 646 3.249 647 1 10 2.0 4.388 505 4.389 062 4.389 056 −1 B ng 1.9. Phương pháp Adams - Moulton áp d ng vào bài toán (2.9); giá tr tiên đoán đư c tính b ng (2.7) và các giá tr đư c hi u chính b ng (2.8).
14 Chương 2 PHƯƠNG PHÁP NHI U NÚT Phương pháp m t bư c(m t nút) là phương pháp mà trong m i bư c chúng ta dùng các k t qu đ t đư c t các bư c trư c m t cách có th t . T t c các phương pháp trong chương 1 đ u là phương pháp m t bư c. Trái l i, phương pháp nhi u nút là phương pháp mà m i bư c dùng các giá tr t các bư c th t . Lý do mà chúng ta dùng các bư c tương đương là thêm các thông tin có th làm tăng đ chính xác c a k t qu . Các phương pháp đó s l n đư c trình bày như dư i đây: 2.1 Phương pháp Adams - Bashforth Ta hãy xét phương trình vi phân có đi u ki n đ u sau y = f (x, y ), y (x0 ) = y0 (2.1) như trư c, v i f như là m t bài toán có nghi m duy nh t trên đo n ch a đi m x0 . Ta l y tích phân y = f (x, y ) t xn đ n xn+1 = xn + h: xn+1 xn+1 y (x)dx = y (xn+1 ) − y (xn ) = f (x, y (x))dx. xn xn Ta thay f (x, y (x)) b i đa th c n i suy (interpolation polynomical) p(x), đ ta đi tính tích phân, đi u này cho k t qu x p xĩ yn+1 c a y (xn+1 ) và yn c a y (xn ) xn+1 yn+1 = yn + p(x)dx (2.2) xn Ch n các giá tr khác nhau c a p(x) s t o ra các phương pháp khác nhau. Ta gi i thích các nguyên lý b ng cách dùng đa th c b c ba, kí hi u p3 (x) t i các đi m cách đ u xn , xn−1 , xn−2 , xn−3 có giá tr l n lư t fn = f (xn , yn ), fn−1 = f (xn−1 , yn−1 ),
15 fn−2 = f (xn−2 , yn−2 ), fn−3 = f (xn−3 , yn−3 ). (2.3) Đi u này s đưa ra các bi u th c t t. Ta có th tính đư c p3 (x) t bi u th c Newton. 1 1 2 3 p3 (x) = fn + r fn + r(r + 1) fn + r(r + 1)(r + 2) fn 2 6 v i r = (x − xn )/h. Ta l y tích phân p3 (x) theo x t xn đ n xn+1 = xn + h, vì v y theo r t 0 đ n 1. Vì x = xn + hr nên dx = hdr. Tích phân c a 1 r(r + 1) là 12 và 5 2 1 3 c a 6 r(r + 1)(r + 2) là 8 . Vì v y ta đ t đư c: xn+1 1 1 5 3 2 3 p3 dx = h P3 dr = h fn + fn + fn + fn (2.4) 2 12 8 xn 0 Ta bi t r ng fn = fn − fn−1 2 fn = fn − 2fn−1 + fn−2 3 fn = fn − 3fn−1 + 3fn−2 − fn−3 . Th vào (2.4) và gom các s h ng l i ta thư đư c bi u th c nhi u bư c c a Phương pháp Adams - Bashforth b c b n: h yn+1 = yn + (55fn − 59fn−1 + 37fn−2 − 9fn−3 ) (2.5) 24 Bi u th c này di n t giá tr m i yn+1 [x p xĩ nghi m y c a (2.1) t i xn+1 ] trong các s c a b n giá tr c a f đư c tính t giá tr c a y đ t đư c trong th t b n bư c. Sai s b ph n là s h5 , có th đư c bi u di n đ sai s toàn ph n là s h4 ; vì v y (2.5) đư c đ nh nghĩa là phương pháp b c b n. 2.2 Các phương pháp Adams - Moulton Các phương pháp này đ t đư c n u cho p(x) trong (2.2) ch n đa th c n i suy f (x, y (x)) t i xn+1 , xn , xn−1 , · · · (trái v i xn , xn−1 , · · · đã đư c dùng trư c đó, đây là các đi m chính). Ta gi i thích quy lu t cho đa th c b c ba p3 (x) n i suy t i các đi m xn+1 , xn , xn−1 , xn−2 . (Trư c đây ta đã có t i các đi m xn , xn−1 , xn−2 , xn−3 . Ta tính đa th c n i suy nhưng lúc này ta thay r = (x − xn+1 h), ta đư c 1 1 2 3 p3 (x) = fn+1 + r fn+1 + r(r + 1) fn+1 + r(r + 1)(r + 2) fn+1 . 2 6 Ta l y tích phân theo x t xn đ n xn+1 như trư c. Như v y, r tương ng t −1 đ n 0. Ta đư c xn+1 1 1 1 2 3 p3 (x)dx = h fn+1 − fn+1 − fn+1 − fn+1 . 2 12 24 xn
16 Ta cũng có bi u th c xn+1 h yn+1 = yn + p3 (x)dx = yn + (9fn+1 + 19fn − 5fn−1 + fn−2 ) (2.6) 24 xn Bi u th c trên thư ng đư c g i là bi u th c Adams - Moulton. Nó là m t đa th c n vì fn+1 = f (xn+1 , yn+1 ) xu t hi n v ph i, vì v y nó đ nh nghĩa yn+1 ch là n, khác v i (2.5), đó là bi u th c n không báo hàm yn+1 v ph i. Đ dùng ∗ (2.6) ta tiên đoán m t giá tr yn+1 , ch ng h n dùng (2.5) là h ∗ yn+1 = yn + (55fn − 59fn−1 + 37fn−2 − 9fn−3 ) (2.7) 24 ∗ ∗ Giá tr đúng m i yn+1 đ t đư c t (2.6) v i fn+1 đư c thay b ng fn+1 = f (xn+1 , yn+1 ) vì v y đ t đư c giá tr như sau: h (9f ∗ + 19fn − 5fn−1 + fn−2 ) yn+1 = yn + (2.8) 24 n+1 Phương pháp predictor - corrector (2.7), (2.8) này đư c g i là phương pháp Adams - Moulton b c b n Phương pháp này ta s d ng qua hao bư c: Bư c (2.7) đư c g i là predictor; Hi u ch nh các giá tr yn+1 trong (2.8) v a tìm đư c b ng cách s d ng các bi u th c n, bư c này g i là corretor và có th l p l i nhi u l n cho đ n khi k t qu đ t đ chính xác như ý. Ví d 2.1. Dùng tiên đoán Adams - Bashforth (2.7), và hi u ch nh Adams - Moultor (2.8) đ gi i phương trình y = x + y, y (0) = 0 (2.9) b ng (2.7) và (2.8) trên đo n 0 ≤ x ≤ 2, ch n h = 0.2 Gi i: Đây là bài toán gi ng như bài toán trong ví d 1.1-1.2 trong chương 1 nên ta có th so sánh các k t qu . Ta tính các giá tr ban đ u y1 , y2 , y3 b ng phương pháp Runge - Kutta c đi n. R i trong m i bư c ta tiên đoán b ng cách dùng (2.7) và hi u ch nh b ng cách dùng (2.8) trư c khi ti n hành các bư c ti p theo. Các k t qu trong b ng 1.9 đư c bi u di n và so sánh v i k t qu chính xác. Ta th y r ng tính chính xác đư c c i ti n, kh o sát chính xác. Nh n xét v s so sánh c a các phương pháp . Bi u th c Adams - Moulton là chính xác hơn bi u th c Adams - Bashforth cùng b c. Phương pháp (2.7); (2.8) là phương pháp s n đ nh, trái l i ngo i tr s d ng (2.7) có th gây ra s b t n đ nh. S đi u ch nh m i bư c s tương đ i đơn gi n. N u |Corrector − P redictor| > T OL, dùng phép n i suy đ phát sinh k t qu cũ t i m t n a t i bư c hi n th i và th h/2 như bư c m i. Phương pháp (2.7), (2.8) ch c n 2 giá tr cho m i bư c, phương pháp Runge - Kutta c n 4 giá tr ; tuy nhiên phương pháp Runge - Kutta có th có s bư c l n hơn hai r t l n.
17 Chương 3 Phương pháp cho h phương trình và các phương trình b c cao Ta xét h phương trình vi phân c p m t có đi u ki n đ u sau: y = f (x, y ), y (x0 ) = y0 , (3.1) như v y y1 = f1 (x, y1 , · · · , ym ) y2 = f2 (x, y1 , · · · , ym ) .... .. ........................ ym = fm (x, y1 , · · · , ym ) Gi s f là bài toán có duy nh t nghi m theo x trên đo n ch a x0 . H phương trình (3.1) g m các phương trình vi phân c p m y (m) = f (x, y, y , y , · · · , y (m−1) ) (3.2) và đi u ki n đ u tương ng là y (x0 ) = K1 , y (x0 ) = K2 , · · · , y (m−1) = Km Ta có th đ t: y3 = y , · · · , ym = y (m−1) y1 = y, y2 = y , (3.3) Thì ta đ t đư c y1 = y2 , y 2 = y3 , · · · , ym−1 = ym , ym = f (x, y1 , · · · , ym ). (3.4) và đi u ki n đ u y1 (x0 ) = K1 , y2 (x0 ) = K2 , · · · , ym (x0 ) = Km . 3.1 Phương pháp Euler cho h phương trình Các phương pháp cho phương trình vi phân b c nh t đơn gi n có th đư c m r ng cho h phương trình (3.1) b ng cách vi t dư i d ng vectơ hàm y và f thay cho các hàm y và f trư c đây, đ y x là m t đ i lư ng vô hư ng.
18 Ví d 3.1. Phương pháp Euler cho phương trình vi ph n c p hai. H kh i lư ng - lò xo. Gi i bài toán t t d n cho h kh i lư ng - lò xo y + 2y + 0.75y = 0, y (0) = 3, y (0) = −2.5 b ng cách dùng phương pháp Euler cho h th ng v i bư c nhãy h = 0.2 và x t 0 đ n 1. Gi i: Phương pháp Euler (1.3) chương 1 khái quát hoá cho h th ng yn+1 = yn + hf (xn , yn ). (3.5) v i các thành ph n y1,n+1 = y1,n + hf1 (xn , y1,n , y2,n ) y2,n+1 = y2,n + hf2 (xn , y1,n , y2,n ) và tương t cho h phương trình nhi u hơn hai phương trình. D a vào (3.4) các phương trình đã cho chuy n đ i thành h y1 = f1 (x, y1 , y2 ) = y2 y2 = f2 (x, y1 , y2 ) = −2y2 − 0.75y1 Do đó (3.5) tr thành y1,n+1 = y1,n + 0.2y2,n y2,n+1 = y2,n + 0.2(−2y2,n − 0.75y1,n ) Các đi u ki n đ u là y (0) = y1 (0) = 3, y (0) = y2 (0) = −2.5 các tính toán đư c trình bày trong b ng 3.1. Như các phương pháp đơn giãn, các k t qu s không đ chính xác cho m c đích tính toán. T t nhiên ví d ch là đ minh ho cho phương pháp còn bài toán này ta có th gi i đư c nghi m m t cách chính xác y = y1 = 2e−0.5x + e−1.5x , y = y2 = −e−0.5x − 1.5e−1.5x y1 Sai s y2 Sai s n xn y1,n chính xác = y1 − y1,n y2,n chính xác = y2 − y2,n 1 2 0 0.0 3.00000 3.00000 0.00000 -2.50000 -2.50000 0.00000 1 0.2 2.50000 2.55049 0.05049 -1.95000 -2.01606 -0.06606 2 0.4 2.11000 2.18627 0.76270 -1.54500 -1.64195 -0.09695 3 0.6 1.80100 1.88821 0.08721 -1.24350 -1.35067 -0.10717 4 0.8 1.55230 1.64183 0.08953 -1.01625 -1.12211 -0.10586 5 1.0 1.34905 1.43619 0.08714 -0.84260 -0.94123 -0.09863 B ng 3.1. Phương pháp Euler cho ví d 3.1
19 3.2 Các phương pháp Runge - Kutta cho h phương trình Đư c dùng cho h phương trình (3.1) b ng cách vi t dư i d ng vectơ, thành ph n m cho m = 1 ép bu c các bi u th c vô hư ng trư c. T phương pháp Runge - Kutta c đi n b c b n trong b ng 1.4 ta đư c y (x0 ) = y0 , giá tr đ u (3.6) và cho n = 0, 1, · · · , N − 1 (N s các bư c), ta đ t đư c các tính ch t ph tr sau 1 1 k1 = hf (xn , yn ), k2 = hf (xn + h, yn + k1 ), 2 2 1 1 k3 = hf (xn + h, yn + k2 ), k4 = hf (xn + h, yn + k3 ) (3.7) 2 2 và giá tr m i [x p xĩ c a nghi m y (x) t i xn+1 = x0 + (n + 1)h] 1 yn+1 = yn + (k1 + 2k2 + 2k3 + k4 ) (3.8) 6 Ví d 3.2. Phương pháp Runge - Kutta cho h phương trình. Phương trình Airy. Hàm Airy Ai(x) Gi i bài toán sau 1 1 y = xy, y (0) = 2/3 = 0.35502805, y (0) = − 1/3 = −0.25881940 3 Γ(2/3) 3 Γ(1/3) b ng phương pháp Runge - Kutta cho h phương trình v i h = 0.2, làm 5 bư c. Đây là phương trình Airy, đư c phát sinh trong quang h c. Γ là hàm gamma. Đi u ki n đ u đ t đư c trong đi u ki n chu n, các hàm Ai(x) đã đư c li t kê và kh o sát. Gi i: Cho y = xy , đ t y1 = y, y2 = y1 = y , ta đư c t h (3.4) y1 = y2 , y2 = xy1 . f2 ]T trong (3.1) có các thành ph n đây f = [f1 f1 (x, y ) = y2 , f2 (x, y ) = xy1 Ta vi t các thành ph n (??). Đi u ki n đ u (3.6) là y1,0 = 0.35508805, y2,0 = −0.25881940. Ta có th vi t k1 = a, k2 = b, k3 = c, k4 = d a2 ]T r i (3.7) có d ng đ a = [a1 y2,n a=h xn y1,n 1 y2,n + 2 a2 b=h (xn + 1 h)(y1,n + 1 a1 ) 2 2
20 y2,n + 1 b2 c = h (3.9) 2 (xn + 2 h)(y1,n + 1 b1 ) 1 2 y2,n + c2 d = h (xn + h)(y1,n + c1 ) Ch ng h n, thành ph n th hai c a b đ t đư c như theo sau, f (x, y ) có thành ph n th hai f2 (x, y ) = xy1 . Trong b(= k2 ) s th nh t là x = xn + 1 h. S th hai 2 trong b là y = yn + 2 a và thành ph n đ u tiên là y1 = y1,n + 1 a1 . Cùng v i, 1 2 xy1 = (xn + 1 h)(y1,n + 1 a1 ). Tương t tính đư c các thành ph n khác trong (3.9). 2 2 Cu i cùng 1 yn+1 = yn + (a + 2b + 2c + d) (3.10) 6 B ng 3.2 di n t các giá tr y (x) = y1 (x) c a hàm Airy Ai(x) và đ o hàm c a nó là y (x) = y2 (x) cũng như sai s c a y (x). 108 × sai s c a y1 n xn y1,n (xn ) y1 (xn ) chính xác y2,n (xn ) 0 0.0 0.35502 805 0.35502 805 0 -0.25881 940 1 0.2 0.30370 303 0.30370 315 12 -0.25240 464 2 0.4 0.25474 211 0.25474 235 24 -0.23583 073 3 0.6 0.20979 973 0.20980 006 33 -0.21279 185 4 0.8 0.16984 596 0.16984 632 36 -0.18641 171 5 1.0 0.13529 207 0.13529 242 35 -0.15914 687 B ng 3.2. Phương pháp Runge - Kutta cho h phương trình: Các giá tr y1,n (xn ) c a hàm Airy Ai(x) trong ví d 3.2 3.3 Phương pháp Runge - Kutta - Nystr¨m (phương o pháp RKN) Phương pháp RKN là m r ng tr c ti p phương pháp RK(phương pháp Runge - Kutta) cho phương trình vi phân c p hai y = f (x, y, y ) như đã cho b ng phép toán E. J. Nystr¨m. Đa th c t t nh t đư c bi t, đây n = 0, 1, · · · , N − 1 (N là o 1 1 s bư c nh y). Các bi u th c dư i v i K = 2 h(yn + 2 k1 ), L = h(yn + k3 ) 1 k1 = hf (xn , yn , yn ) 2 1 1 k2 = hf (xn + h, yn + K, yn + k1 ) 2 2 1 1 k3 = hf (xn + h, yn + K, yn + k2 ) (3.11) 2 2 1 k4 = hf (xn + h, yn + L, yn + 2k3 ) 2