›
›

Bài giảng Lập trình Matlab - Chương 4: Ứng dụng Matlab trong Giải tích số

pdf 104 trang Gia Huy 17/05/2022 3350
Download
Bạn đang xem 20 trang mẫu của tài liệu "Bài giảng Lập trình Matlab - Chương 4: Ứng dụng Matlab trong Giải tích số", để tải tài liệu gốc về máy bạn click vào nút DOWNLOAD ở trên

Tài liệu đính kèm:

  • pdfbai_giang_lap_trinh_matlab_chuong_4_ung_dung_matlab_trong_gi.pdf

Nội dung text: Bài giảng Lập trình Matlab - Chương 4: Ứng dụng Matlab trong Giải tích số

  1. Matlab trong Giải tích số 1/57 Chương 4: Ứng dụng Matlab trong Giải tích số Viện Toán ứng dụng và Tin học, ĐHBK Hà Nội Hà Nội, tháng 8 năm 2015
  2. Matlab trong Giải tích số 2/57 Đa thức nội suy Nội dung 1 Đa thức nội suy Nội suy Lagrange Nội suy Newton Nội suy bằng đa thức Chebyshev Nội suy bằng đa thức Hermit (đọc thêm) 2 Giải gần đúng phương trình Phương pháp chia đôi Phương pháp dây cung Phương pháp Newton - Raphson 3 Giải gần đúng hệ phương trình Phương pháp lặp đơn Phương pháp Jacobi Phương pháp Gauss-Seidel Phương pháp Newton giải hệ phương trình phi tuyến 4 Giải gần đúng phương trình vi phân thường Phương pháp Euler Phương pháp Euler cải tiến (Modified Euler hay RK-2) Phương pháp Runge-Kutta Hệ phương trình vi phân thường và phương trình vi phân cấp cao
  3. Matlab trong Giải tích số 3/57 Đa thức nội suy Đa thức nội suy Trong thực tế, nhiều khi ta phải tìm hàm = ( ) mà chỉ biết giá trị 푖 tại các điểm 푖 ∈ [ , ](푖 = 0, 1, . . . , 푛). Hoặc trong nhiều trường hợp biểu diễn giải tích của ( ) đã cho nhưng quá cồng kềnh. Khi dùng phép nội suy ta có thể dễ dàng tính được tại bất kỳ điểm ∈ [ , ] mà độ chính xác không kém bao nhiêu. Bài toán đặt ra: Cho các mốc nội suy ≤ 0 < 1 < ··· < 푛 ≤ . Hãy tìm đa thức (bậc 푛) 푛 ∑︀ 푖 푃푛( ) = 푖 sao cho: 푖=0 푃푛( 푖) = 푖 = ( 푖)(푖 = 0 ÷ 푛) (1.1) Đa thức 푃푛( ) gọi là đa thức nội suy của hàm = ( ). Ta chọn đa thức để nội suy hàm vì đa thức là loại hàm đơn giản, luôn có đạo hàm và nguyên hàm. Việc tính giá trị của nó theo thuật toán Horner cũng đơn giản.
  4. Matlab trong Giải tích số 4/57 Đa thức nội suy Đa thức nội suy Cách tiếp cận Vandermond Các hệ số 0, 1, . . . , 푛 của đa thức nội suy bậc 푛 có thể được tính bằng cách giải hệ ⎡ 푛 푛−1 2 ⎤ ⎡ ⎤ ⎡ ⎤ 0 0 ··· 0 0 1 0 0 푛 푛−1 2 ⎢ 1 1 ··· 1 1 1⎥ ⎢ 1 ⎥ ⎢ 1 ⎥ ⎢ 푛 푛−1 2 ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ 2 2 ··· 2 2 1⎥ ⎢ 2 ⎥ ⎢ 2 ⎥ ⎢ ⎥ ⎢ ⎥ = ⎢ ⎥ ⎢ . . ⎥ ⎢ . ⎥ ⎢ . ⎥ ⎣ . . ⎦ ⎣ . ⎦ ⎣ . ⎦ 푛 푛−1 2 푛 푛 ··· 푛 푛 1 푛 푛 hay = 푛 ∏︀ Hệ trên có định thức Vandermond | | = ( 푗 − 푖) ̸= 0 nên có nghiệm 1≤푖<푗 duy nhất.
  5. Matlab trong Giải tích số 5/57 Đa thức nội suy Nội suy Lagrange Nội suy Lagrange Trước hết tìm đa thức 퐿푖( ) có bậc 푛 sao cho: {︃ 1 nếu 푖 = 푗 퐿푖( 푗 ) = , (∀푖, 푗 = 0 ÷ 푛) 0 nếu 푖 ̸= 푗 Dễ thấy 퐿푖( ) có dạng: ∏︀ ( − 푗 ) 푗̸=푖 퐿푖( ) = ∏︀ ( 푖 − 푗 ) 푗̸=푖 ∏︀ ( − 푗 ) 푛 푛 ∑︀ ∑︀ 푗̸=푖 Đặt 푃 ( ) = 푖.퐿푖( ) = 푖 ∏︀ ta có ngay: 푖=0 푖=1 ( 푖 − 푗 ) 푗̸=푖 푛 ∑︁ 푃 ( 푗 ) = 푖.퐿푖( 푗 ) = 푗 (푗 = 0 ÷ 푛) (1.2) 푖=1 Vậy 푃 ( ) là đa thức nội suy duy nhất cần tìm.
  6. Matlab trong Giải tích số 6/57 Đa thức nội suy Nội suy Newton Nội suy Newton Nội suy Newton tiến Công thức nội suy Lagrange có ưu điểm là đơn giản, dễ lập trình nhưng nếu thêm mốc nội suy thì phải tính lại toàn bộ. Nhược điểm này sẽ được khắc phục trong công thức Newton. Công thức nội suy Newton tiến 푃푛( ) = ( 0) + ( − 0) [ 0, 1] + ( − 0)( − 1) [ 0, 1, 2] + ··· + ( − 0)( − 1) ( − 푛−1) [ 0, 1, . . . , 푛] , trong đó các tỷ hiệu được tính theo công thức ( 푖) − ( 푖−1) [ 푖−1, 푖] = ; 푖 − 푖−1 [ 1, . . . , 푛] − [ 0, . . . , 푛−1] [ 0, 1, . . . , 푛] = (1.3) 푛 − 0
  7. Matlab trong Giải tích số 6/57 Đa thức nội suy Nội suy Newton Nội suy Newton Nội suy Newton tiến Công thức nội suy Lagrange có ưu điểm là đơn giản, dễ lập trình nhưng nếu thêm mốc nội suy thì phải tính lại toàn bộ. Nhược điểm này sẽ được khắc phục trong công thức Newton. Công thức nội suy Newton tiến 푃푛( ) = ( 0) + ( − 0) [ 0, 1] + ( − 0)( − 1) [ 0, 1, 2] + ··· + ( − 0)( − 1) ( − 푛−1) [ 0, 1, . . . , 푛] , trong đó các tỷ hiệu được tính theo công thức ( 푖) − ( 푖−1) [ 푖−1, 푖] = ; 푖 − 푖−1 [ 1, . . . , 푛] − [ 0, . . . , 푛−1] [ 0, 1, . . . , 푛] = (1.3) 푛 − 0
  8. Matlab trong Giải tích số 7/57 Đa thức nội suy Nội suy Newton Nội suy Newton Nội suy Newton lùi Nếu các mốc nội suy được sắp xếp theo thứ tự giảm dần 푛, 푛−1, . . . , 1, 0 thì ta có công thức nội suy Newton lùi xuất phát từ mốc 푛: 푃푛( ) = ( 푛) + ( − 푛) [ 푛, 푛−1] + ( − 푛)( − 푛−1) [ 푛, 푛−1, 푛−2] + ··· + ( − 푛)( − 푛−1) ( − 1) [ 푛, 푛−1, . . . , 1, 0] , trong đó các tỷ hiệu được tính như trong công thức (1.3).
  9. Matlab trong Giải tích số 8/57 Đa thức nội suy Nội suy Newton Sai số của phép nội suy Định lý 1.1 Giả sử hàm : R → R khả vi liên tục đến cấp 푛 + 1 trên [ , ] (푛+1) ( ∈ [ , ]) và 푖 ∈ [ , ], 푖 = 0 : 푛. Khi đó tồn tại 휉 = 휉( ) ∈ [ , ] sao cho 1 ( ) − 푃 ( ) = (푛+1)(휉)( − ) ( − ). 푛 (푛 + 1)! 0 푛 Từ đó ta có công thức ước lượng sai số 1 | ( ) − 푃 ( )| ≤ ( ) |( − ) ( − )| , 푛 (푛 + 1)! 푛+1 0 푛 ⃒ (푛+1) ⃒ trong đó 푛+1( ) = max ⃒ ( )⃒. ∈[ , ] ⃒ ⃒
  10. Matlab trong Giải tích số 9/57 Đa thức nội suy Nội suy bằng đa thức Chebyshev Nội suy bằng đa thức Chebyshev Xét trường hợp nội suy đa thức cho hàm ( ) trên đoạn [−1, 1] dựa trên các mốc nội suy −1 ≤ 0 < 1 < ··· < 푛 ≤ 1. Khi đó công thức đánh giá sai số của các đa thức nội suy Lagrange và Newton đều có dạng 1 ⃒ (푛+1) ⃒ | ( ) − 푃푛( )| ≤ |푤( )| max ⃒ ( )⃒ , (푛 + 1)! ∈[−1,1] ⃒ ⃒ trong đó đa thức bậc 푛 + 1: 푤( ) = ( − 0)( − 1) ( − 푛). 푛 Ta muốn chọn các mốc nội suy { 푖} để cực tiểu giá trị max |푤( )|. 푖=0 −1≤ ≤1 Điều này dẫn tới việc sử dụng đa thức nội suy Chebyshev.
  11. Matlab trong Giải tích số 9/57 Đa thức nội suy Nội suy bằng đa thức Chebyshev Nội suy bằng đa thức Chebyshev Xét trường hợp nội suy đa thức cho hàm ( ) trên đoạn [−1, 1] dựa trên các mốc nội suy −1 ≤ 0 < 1 < ··· < 푛 ≤ 1. Khi đó công thức đánh giá sai số của các đa thức nội suy Lagrange và Newton đều có dạng 1 ⃒ (푛+1) ⃒ | ( ) − 푃푛( )| ≤ |푤( )| max ⃒ ( )⃒ , (푛 + 1)! ∈[−1,1] ⃒ ⃒ trong đó đa thức bậc 푛 + 1: 푤( ) = ( − 0)( − 1) ( − 푛). 푛 Ta muốn chọn các mốc nội suy { 푖} để cực tiểu giá trị max |푤( )|. 푖=0 −1≤ ≤1 Điều này dẫn tới việc sử dụng đa thức nội suy Chebyshev.
  12. Matlab trong Giải tích số 10/57 Đa thức nội suy Nội suy bằng đa thức Chebyshev Nội suy bằng đa thức Chebyshev Định nghĩa 1.1 Các hàm 푛( ) = cos(푛 arccos( )), 푛 = 0, 1, 2, gọi là các đa thức Chebyshev trong đoạn [−1, 1]. Chú ý 1.1 Các hàm trên thực sự là các đa thức. Thật vậy, đặt 훿 = arccos( ). Đồng nhất thức cos(푛 + 1)훿 + cos(푛 − 1)훿 = 2 cos 훿 cos 푛훿 cho ta công thức truy hồi 푛+1( ) = 2 푛( ) − 푛−1( ). Với 0( ) = 1, 1( ) = , rõ ràng 푛( ) ∈ 풫푛.
  13. Matlab trong Giải tích số 11/57 Đa thức nội suy Nội suy bằng đa thức Chebyshev Nội suy bằng đa thức Chebyshev Bảng một số đa thức nội suy Chebyshev đầu tiên 0( ) = 1 1( ) = 2 2( ) = 2 − 1 3 3( ) = 4 − 3 4 2 4( ) = 8 − 8 − 1 5 3 5( ) = 16 − 20 + 5 6 4 2 6( ) = 32 − 48 + 18 − 1 7 5 3 7( ) = 64 − 112 + 56 − 7
  14. Matlab trong Giải tích số 12/57 Đa thức nội suy Nội suy bằng đa thức Chebyshev Nội suy bằng đa thức Chebyshev Một số tính chất của đa thức Chebyshev Các hệ số của đa thức Chebyshev 푛 đều nguyên. 푛−1 Hệ số ứng với bậc cao nhất là 푛 = 2 . 2푛 là hàm chẵn, 2푛+1 là hàm lẻ. | 푛( )| ≤ 1 với ∈ [−1, 1] và 푛( ) = 1 với = cos( /푛). 푛 푛(1) = 1, 푛(−1) = (−1) . (︂ (2 − 1) )︂ ( ) = 0 với = cos , = 1, . . . , 푛. 푛 2
  15. Matlab trong Giải tích số 12/57 Đa thức nội suy Nội suy bằng đa thức Chebyshev Nội suy bằng đa thức Chebyshev Một số tính chất của đa thức Chebyshev Các hệ số của đa thức Chebyshev 푛 đều nguyên. 푛−1 Hệ số ứng với bậc cao nhất là 푛 = 2 . 2푛 là hàm chẵn, 2푛+1 là hàm lẻ. | 푛( )| ≤ 1 với ∈ [−1, 1] và 푛( ) = 1 với = cos( /푛). 푛 푛(1) = 1, 푛(−1) = (−1) . (︂ (2 − 1) )︂ ( ) = 0 với = cos , = 1, . . . , 푛. 푛 2
  16. Matlab trong Giải tích số 12/57 Đa thức nội suy Nội suy bằng đa thức Chebyshev Nội suy bằng đa thức Chebyshev Một số tính chất của đa thức Chebyshev Các hệ số của đa thức Chebyshev 푛 đều nguyên. 푛−1 Hệ số ứng với bậc cao nhất là 푛 = 2 . 2푛 là hàm chẵn, 2푛+1 là hàm lẻ. | 푛( )| ≤ 1 với ∈ [−1, 1] và 푛( ) = 1 với = cos( /푛). 푛 푛(1) = 1, 푛(−1) = (−1) . (︂ (2 − 1) )︂ ( ) = 0 với = cos , = 1, . . . , 푛. 푛 2
  17. Matlab trong Giải tích số 12/57 Đa thức nội suy Nội suy bằng đa thức Chebyshev Nội suy bằng đa thức Chebyshev Một số tính chất của đa thức Chebyshev Các hệ số của đa thức Chebyshev 푛 đều nguyên. 푛−1 Hệ số ứng với bậc cao nhất là 푛 = 2 . 2푛 là hàm chẵn, 2푛+1 là hàm lẻ. | 푛( )| ≤ 1 với ∈ [−1, 1] và 푛( ) = 1 với = cos( /푛). 푛 푛(1) = 1, 푛(−1) = (−1) . (︂ (2 − 1) )︂ ( ) = 0 với = cos , = 1, . . . , 푛. 푛 2
  18. Matlab trong Giải tích số 12/57 Đa thức nội suy Nội suy bằng đa thức Chebyshev Nội suy bằng đa thức Chebyshev Một số tính chất của đa thức Chebyshev Các hệ số của đa thức Chebyshev 푛 đều nguyên. 푛−1 Hệ số ứng với bậc cao nhất là 푛 = 2 . 2푛 là hàm chẵn, 2푛+1 là hàm lẻ. | 푛( )| ≤ 1 với ∈ [−1, 1] và 푛( ) = 1 với = cos( /푛). 푛 푛(1) = 1, 푛(−1) = (−1) . (︂ (2 − 1) )︂ ( ) = 0 với = cos , = 1, . . . , 푛. 푛 2
  19. Matlab trong Giải tích số 12/57 Đa thức nội suy Nội suy bằng đa thức Chebyshev Nội suy bằng đa thức Chebyshev Một số tính chất của đa thức Chebyshev Các hệ số của đa thức Chebyshev 푛 đều nguyên. 푛−1 Hệ số ứng với bậc cao nhất là 푛 = 2 . 2푛 là hàm chẵn, 2푛+1 là hàm lẻ. | 푛( )| ≤ 1 với ∈ [−1, 1] và 푛( ) = 1 với = cos( /푛). 푛 푛(1) = 1, 푛(−1) = (−1) . (︂ (2 − 1) )︂ ( ) = 0 với = cos , = 1, . . . , 푛. 푛 2
  20. Matlab trong Giải tích số 13/57 Đa thức nội suy Nội suy bằng đa thức Chebyshev Nội suy bằng đa thức Chebyshev Định lý 1.2 Giả sử 푛 cố định. Trong số tất cả các cách chọn 푤( ) và các mốc phân biệt 푛 푛 { 푖}푖=0 ∈ [−1, 1], đa thức ( ) = 푛+1( )/2 là sự lựa chọn suy nhất thỏa mãn max {| ( )|} ≤ max {|푤( )|} . −1≤ ≤1 −1≤ ≤1 Hơn nữa 1 max {| ( )|} = . −1≤ ≤1 2푛
  21. Matlab trong Giải tích số 14/57 Đa thức nội suy Nội suy bằng đa thức Chebyshev Nội suy bằng đa thức Chebyshev Các mốc nội suy trên đoạn [−1, 1] được xác định bởi (︂ (2푛 + 1 − 2 ) )︂ 푡 = cos , = 0, 1, . . . , 푛 2푛 + 2 Sử dụng phép đổi biến [−1, 1] → [ , ]: (︂ − )︂ + − = 푡 + ⇐⇒ 푡 = 2 − 1 2 2 − Các mốc nội suy trên đoạn [ , ] bất kỳ (︂ − )︂ + = 푡 + 2 2 (︂ (2푛 + 1 − 2 ) )︂ (︂ − )︂ + = cos + , = 0, 1, . . . , 푛. (1.4) 2푛 + 2 2 2
  22. Matlab trong Giải tích số 14/57 Đa thức nội suy Nội suy bằng đa thức Chebyshev Nội suy bằng đa thức Chebyshev Các mốc nội suy trên đoạn [−1, 1] được xác định bởi (︂ (2푛 + 1 − 2 ) )︂ 푡 = cos , = 0, 1, . . . , 푛 2푛 + 2 Sử dụng phép đổi biến [−1, 1] → [ , ]: (︂ − )︂ + − = 푡 + ⇐⇒ 푡 = 2 − 1 2 2 − Các mốc nội suy trên đoạn [ , ] bất kỳ (︂ − )︂ + = 푡 + 2 2 (︂ (2푛 + 1 − 2 ) )︂ (︂ − )︂ + = cos + , = 0, 1, . . . , 푛. (1.4) 2푛 + 2 2 2
  23. Matlab trong Giải tích số 14/57 Đa thức nội suy Nội suy bằng đa thức Chebyshev Nội suy bằng đa thức Chebyshev Các mốc nội suy trên đoạn [−1, 1] được xác định bởi (︂ (2푛 + 1 − 2 ) )︂ 푡 = cos , = 0, 1, . . . , 푛 2푛 + 2 Sử dụng phép đổi biến [−1, 1] → [ , ]: (︂ − )︂ + − = 푡 + ⇐⇒ 푡 = 2 − 1 2 2 − Các mốc nội suy trên đoạn [ , ] bất kỳ (︂ − )︂ + = 푡 + 2 2 (︂ (2푛 + 1 − 2 ) )︂ (︂ − )︂ + = cos + , = 0, 1, . . . , 푛. (1.4) 2푛 + 2 2 2
  24. Matlab trong Giải tích số 15/57 Đa thức nội suy Nội suy bằng đa thức Chebyshev Nội suy bằng đa thức Chebyshev Định lý 1.3 Giả sử 푃푛( ) là đa thức nội suy Lagrange với các mốc nội suy (1.4). Khi đó nếu ∈ 푛+1 [ , ] thì ta có 푛+1 2( − ) {︁⃒ (푛+1) ⃒}︁ | ( ) − 푃푛( )| ≤ max ⃒ ( )⃒ . 4푛+1(푛 + 1)! ≤ ≤ ⃒ ⃒ Ví dụ 1 [︁ ]︁ Xét hàm ( ) = sin trên 0, . Các mốc nội suy Chebyshev 4 (︂ (11 − 2 ) )︂ = cos + , = 0, 1, , 5. 12 8 8 ⃒ ⃒ ⃒ (6) ⃒ −1/2 Sử dụng đánh giá ⃒ ( )⃒ ≤ |− sin ( /4)| = 2 =: ta thu được 6 (︂ )︂ (︁ )︁ 2 −1/2 | ( ) − 푃5( )| ≤ 2 = 0.00000720. 8 6!
  25. Matlab trong Giải tích số 16/57 Đa thức nội suy Nội suy bằng đa thức Chebyshev Nội suy bằng đa thức Chebyshev Công thức nội suy Chebyshev Hàm ( ) được xấp xỉ bởi 푛 ∑︁ (︀ ′)︀ ( ) ≈ 푗 (1.5) 푗=0 trong đó 2 (︂ + )︂ ′ = − ; − 2 푛 푛 1 ∑︀ ′ 1 ∑︀ 0 = ( ) 0 ( ) = ( ); 푛 + 1 =0 푛 + 1 =0 푛 2 ∑︀ ′ 푗 = ( ) 푗 ( ) = 푛 + 1 =0 푛 2 ∑︀ 푗(2푛 + 1 − 2 ) ( ) cos , 푗 = 1, 2, . . . , 푛. 푛 + 1 =0 2(푛 + 1)
  26. Matlab trong Giải tích số 16/57 Đa thức nội suy Nội suy bằng đa thức Chebyshev Nội suy bằng đa thức Chebyshev Công thức nội suy Chebyshev Hàm ( ) được xấp xỉ bởi 푛 ∑︁ (︀ ′)︀ ( ) ≈ 푗 (1.5) 푗=0 trong đó 2 (︂ + )︂ ′ = − ; − 2 푛 푛 1 ∑︀ ′ 1 ∑︀ 0 = ( ) 0 ( ) = ( ); 푛 + 1 =0 푛 + 1 =0 푛 2 ∑︀ ′ 푗 = ( ) 푗 ( ) = 푛 + 1 =0 푛 2 ∑︀ 푗(2푛 + 1 − 2 ) ( ) cos , 푗 = 1, 2, . . . , 푛. 푛 + 1 =0 2(푛 + 1)
  27. Matlab trong Giải tích số 16/57 Đa thức nội suy Nội suy bằng đa thức Chebyshev Nội suy bằng đa thức Chebyshev Công thức nội suy Chebyshev Hàm ( ) được xấp xỉ bởi 푛 ∑︁ (︀ ′)︀ ( ) ≈ 푗 (1.5) 푗=0 trong đó 2 (︂ + )︂ ′ = − ; − 2 푛 푛 1 ∑︀ ′ 1 ∑︀ 0 = ( ) 0 ( ) = ( ); 푛 + 1 =0 푛 + 1 =0 푛 2 ∑︀ ′ 푗 = ( ) 푗 ( ) = 푛 + 1 =0 푛 2 ∑︀ 푗(2푛 + 1 − 2 ) ( ) cos , 푗 = 1, 2, . . . , 푛. 푛 + 1 =0 2(푛 + 1)
  28. Matlab trong Giải tích số 17/57 Đa thức nội suy Nội suy bằng đa thức Chebyshev Nội suy bằng đa thức Chebyshev Ví dụ 2 Tìm đa thức nội suy Chebyshev bậc 3 của hàm ( ) = 푒 trên [−1, 1]. (︂ (2 + 1) )︂ Các mốc nội suy: = cos , = 0, 1, 2, 3; 8 Các hệ số: 3 3 1 ∑︁ 1 ∑︁ = 푒 ( ) = 푒 = 1.26606568 0 4 0 4 =0 =0 3 3 1 ∑︁ 1 ∑︁ = 푒 ( ) = 푒 = 1.13031500 1 2 1 4 =0 =0 3 3 (︂ )︂ 1 ∑︁ 1 ∑︁ 2 + 1 = 푒 ( ) = 푒 cos 2 = 0.27145036 2 2 2 4 8 =0 =0 3 3 (︂ )︂ 1 ∑︁ 1 ∑︁ 2 + 1 = 푒 ( ) = 푒 cos 3 = 0.04379392. 3 24 2 8 =0 =0
  29. Matlab trong Giải tích số 17/57 Đa thức nội suy Nội suy bằng đa thức Chebyshev Nội suy bằng đa thức Chebyshev Ví dụ 2 Tìm đa thức nội suy Chebyshev bậc 3 của hàm ( ) = 푒 trên [−1, 1]. (︂ (2 + 1) )︂ Các mốc nội suy: = cos , = 0, 1, 2, 3; 8 Các hệ số: 3 3 1 ∑︁ 1 ∑︁ = 푒 ( ) = 푒 = 1.26606568 0 4 0 4 =0 =0 3 3 1 ∑︁ 1 ∑︁ = 푒 ( ) = 푒 = 1.13031500 1 2 1 4 =0 =0 3 3 (︂ )︂ 1 ∑︁ 1 ∑︁ 2 + 1 = 푒 ( ) = 푒 cos 2 = 0.27145036 2 2 2 4 8 =0 =0 3 3 (︂ )︂ 1 ∑︁ 1 ∑︁ 2 + 1 = 푒 ( ) = 푒 cos 3 = 0.04379392. 3 24 2 8 =0 =0
  30. Matlab trong Giải tích số 18/57 Đa thức nội suy Nội suy bằng đa thức Chebyshev Nội suy bằng đa thức Chebyshev Ví dụ 2 (tiếp) Đa thức nội suy Chebyshev bậc 3 3 ∑︁ 푃3( ) = 푗 푗 ( ) 푗=0 = 1.26606568 0( ) + 1.13031500 1( ) + 0.27145036 2( ) + 0.04379392 3( ) = 0.99461532 + 0.99893324 + 0.54290072 2 + 0.17517568 3.
  31. Matlab trong Giải tích số 19/57 Đa thức nội suy Nội suy bằng đa thức Hermit (đọc thêm) Nội suy bằng đa thức Hermit Trong một số trường hợp, ta cần tìm hàm đa thức không những đi qua những điểm cho trước mà còn phải thỏa mãn điều kiện về đạo hàm tại các điểm đó. Ta gọi các đa thức đó là đa thức nội suy Hermit. Để đơn giản, ta khảo sát đa thức bậc 3: 3 2 ℎ( ) = 3 + 2 + 1 + 0 ′ ′ đi qua hai điểm ( 0, 0) , ( 1, 1) và có các đạo hàm là 0, 1. Như vậy ta phải tìm các hệ số 푖, 푖 = 0, 3 bằng cách giải hệ phương trình ⎧ℎ ( ) = 3 + 2 + + = ⎪ 0 3 0 2 0 1 0 0 0 ⎪ 3 2 ⎨ℎ ( 1) = 3 1 + 2 1 + 1 1 + 0 = 1 ′ ′ ℎ ( ) = 3 2 + 2 + = ⎪ 0 3 0 2 0 1 0 ⎩⎪ ′ 2 ′ ℎ ( 1) = 3 3 1 + 2 2 1 + 1 = 1
  32. Matlab trong Giải tích số 20/57 Đa thức nội suy Nội suy bằng đa thức Hermit (đọc thêm) Nội suy bằng đa thức Hermit Các đạo hàm bậc nhất được tính gần đúng bởi ℎ ( + 휀) − ℎ ( ) − ′ = 0 0 =: 2 0 0 휀 휀 ℎ ( ) − ℎ ( − 휀) − ′ = 1 1 =: 1 3 0 휀 휀 Bây giờ ta tìm đa thức nội suy Lagange hay Newton đi qua 4 điểm (︀ ′ )︀ (︀ ′ )︀ ( 0, 0) , 2 = 0 + 휀, 2 = 0 + 0휀 , 3 = 1 − 휀, 3 = 1 + 1휀 , ( 1, 1) .
  33. Matlab trong Giải tích số 21/57 Giải gần đúng phương trình Nội dung 1 Đa thức nội suy Nội suy Lagrange Nội suy Newton Nội suy bằng đa thức Chebyshev Nội suy bằng đa thức Hermit (đọc thêm) 2 Giải gần đúng phương trình Phương pháp chia đôi Phương pháp dây cung Phương pháp Newton - Raphson 3 Giải gần đúng hệ phương trình Phương pháp lặp đơn Phương pháp Jacobi Phương pháp Gauss-Seidel Phương pháp Newton giải hệ phương trình phi tuyến 4 Giải gần đúng phương trình vi phân thường Phương pháp Euler Phương pháp Euler cải tiến (Modified Euler hay RK-2) Phương pháp Runge-Kutta Hệ phương trình vi phân thường và phương trình vi phân cấp cao
  34. Matlab trong Giải tích số 22/57 Giải gần đúng phương trình Giải gần đúng phương trình Trong chương này, chúng ta sẽ nghiên cứu một số phương pháp giải phương trình một biến số: ( ) = 0, (2.1) trong đó ( ) là một hàm số (hay siêu việt). Phương trình (2.1) chỉ giải được đúng trong một số trường hợp đặc biệt, nói chung rất phức tạp, do đó chúng ta phải tìm cách giải gần đúng. Ngoài ra các hệ số của ( ) trong thực tế chỉ biết gần đúng vì thế việc giải đúng (2.1) chẳng những không thực hiện được mà nhiều khi không có ý nghĩa. Để giải (2.1) thông thường có hai bước: 1 Giải sơ bộ: Đi tìm một khoảng đủ bé chứa nghiệm Vây nghiệm: Tìm đoạn bé chứa các nghiệm Tách nghiệm: Tách các đoạn bé, mỗi đoạn chỉ chứa một nghiệm 2 Giải kiện toàn: Tìm nghiệm với độ chính xác cần thiết
  35. Matlab trong Giải tích số 22/57 Giải gần đúng phương trình Giải gần đúng phương trình Trong chương này, chúng ta sẽ nghiên cứu một số phương pháp giải phương trình một biến số: ( ) = 0, (2.1) trong đó ( ) là một hàm số (hay siêu việt). Phương trình (2.1) chỉ giải được đúng trong một số trường hợp đặc biệt, nói chung rất phức tạp, do đó chúng ta phải tìm cách giải gần đúng. Ngoài ra các hệ số của ( ) trong thực tế chỉ biết gần đúng vì thế việc giải đúng (2.1) chẳng những không thực hiện được mà nhiều khi không có ý nghĩa. Để giải (2.1) thông thường có hai bước: 1 Giải sơ bộ: Đi tìm một khoảng đủ bé chứa nghiệm Vây nghiệm: Tìm đoạn bé chứa các nghiệm Tách nghiệm: Tách các đoạn bé, mỗi đoạn chỉ chứa một nghiệm 2 Giải kiện toàn: Tìm nghiệm với độ chính xác cần thiết
  36. Matlab trong Giải tích số 22/57 Giải gần đúng phương trình Giải gần đúng phương trình Trong chương này, chúng ta sẽ nghiên cứu một số phương pháp giải phương trình một biến số: ( ) = 0, (2.1) trong đó ( ) là một hàm số (hay siêu việt). Phương trình (2.1) chỉ giải được đúng trong một số trường hợp đặc biệt, nói chung rất phức tạp, do đó chúng ta phải tìm cách giải gần đúng. Ngoài ra các hệ số của ( ) trong thực tế chỉ biết gần đúng vì thế việc giải đúng (2.1) chẳng những không thực hiện được mà nhiều khi không có ý nghĩa. Để giải (2.1) thông thường có hai bước: 1 Giải sơ bộ: Đi tìm một khoảng đủ bé chứa nghiệm Vây nghiệm: Tìm đoạn bé chứa các nghiệm Tách nghiệm: Tách các đoạn bé, mỗi đoạn chỉ chứa một nghiệm 2 Giải kiện toàn: Tìm nghiệm với độ chính xác cần thiết
  37. Matlab trong Giải tích số 22/57 Giải gần đúng phương trình Giải gần đúng phương trình Trong chương này, chúng ta sẽ nghiên cứu một số phương pháp giải phương trình một biến số: ( ) = 0, (2.1) trong đó ( ) là một hàm số (hay siêu việt). Phương trình (2.1) chỉ giải được đúng trong một số trường hợp đặc biệt, nói chung rất phức tạp, do đó chúng ta phải tìm cách giải gần đúng. Ngoài ra các hệ số của ( ) trong thực tế chỉ biết gần đúng vì thế việc giải đúng (2.1) chẳng những không thực hiện được mà nhiều khi không có ý nghĩa. Để giải (2.1) thông thường có hai bước: 1 Giải sơ bộ: Đi tìm một khoảng đủ bé chứa nghiệm Vây nghiệm: Tìm đoạn bé chứa các nghiệm Tách nghiệm: Tách các đoạn bé, mỗi đoạn chỉ chứa một nghiệm 2 Giải kiện toàn: Tìm nghiệm với độ chính xác cần thiết
  38. Matlab trong Giải tích số 23/57 Giải gần đúng phương trình Phương pháp chia đôi Phương pháp chia đôi Giả sử hàm số ( ) liên tục trên đoạn [ , ] và ( ). ( ) < 0. + Chia [ , ] thành 2 phần bởi điểm giữa = 2 1 Nếu ( ) = 0 thì nghiệm 휉 = 2 Nếu ( ) ̸= 0 thì chọn [ , ] hoặc [ , ] mà giá trị hàm tại hai đầu trái dấu và ký hiệu là [ 1, 1]. Đối với [ 1, 1] lại tiến hành như [ , ]. Ưu điểm của phương pháp chia đôi là thuật toán đơn giản do đó dễ lập trình. Tuy nhiên do phương pháp chia đôi sử dụng rất ít thông tin về hàm ( ) nên tốc độ hội tụ khá chậm và chỉ sử dụng để giải sơ bộ phương trình.
  39. Matlab trong Giải tích số 23/57 Giải gần đúng phương trình Phương pháp chia đôi Phương pháp chia đôi Giả sử hàm số ( ) liên tục trên đoạn [ , ] và ( ). ( ) < 0. + Chia [ , ] thành 2 phần bởi điểm giữa = 2 1 Nếu ( ) = 0 thì nghiệm 휉 = 2 Nếu ( ) ̸= 0 thì chọn [ , ] hoặc [ , ] mà giá trị hàm tại hai đầu trái dấu và ký hiệu là [ 1, 1]. Đối với [ 1, 1] lại tiến hành như [ , ]. Ưu điểm của phương pháp chia đôi là thuật toán đơn giản do đó dễ lập trình. Tuy nhiên do phương pháp chia đôi sử dụng rất ít thông tin về hàm ( ) nên tốc độ hội tụ khá chậm và chỉ sử dụng để giải sơ bộ phương trình.
  40. Matlab trong Giải tích số 23/57 Giải gần đúng phương trình Phương pháp chia đôi Phương pháp chia đôi Giả sử hàm số ( ) liên tục trên đoạn [ , ] và ( ). ( ) < 0. + Chia [ , ] thành 2 phần bởi điểm giữa = 2 1 Nếu ( ) = 0 thì nghiệm 휉 = 2 Nếu ( ) ̸= 0 thì chọn [ , ] hoặc [ , ] mà giá trị hàm tại hai đầu trái dấu và ký hiệu là [ 1, 1]. Đối với [ 1, 1] lại tiến hành như [ , ]. Ưu điểm của phương pháp chia đôi là thuật toán đơn giản do đó dễ lập trình. Tuy nhiên do phương pháp chia đôi sử dụng rất ít thông tin về hàm ( ) nên tốc độ hội tụ khá chậm và chỉ sử dụng để giải sơ bộ phương trình.
  41. Matlab trong Giải tích số 23/57 Giải gần đúng phương trình Phương pháp chia đôi Phương pháp chia đôi Giả sử hàm số ( ) liên tục trên đoạn [ , ] và ( ). ( ) < 0. + Chia [ , ] thành 2 phần bởi điểm giữa = 2 1 Nếu ( ) = 0 thì nghiệm 휉 = 2 Nếu ( ) ̸= 0 thì chọn [ , ] hoặc [ , ] mà giá trị hàm tại hai đầu trái dấu và ký hiệu là [ 1, 1]. Đối với [ 1, 1] lại tiến hành như [ , ]. Ưu điểm của phương pháp chia đôi là thuật toán đơn giản do đó dễ lập trình. Tuy nhiên do phương pháp chia đôi sử dụng rất ít thông tin về hàm ( ) nên tốc độ hội tụ khá chậm và chỉ sử dụng để giải sơ bộ phương trình.
  42. Matlab trong Giải tích số 23/57 Giải gần đúng phương trình Phương pháp chia đôi Phương pháp chia đôi Giả sử hàm số ( ) liên tục trên đoạn [ , ] và ( ). ( ) < 0. + Chia [ , ] thành 2 phần bởi điểm giữa = 2 1 Nếu ( ) = 0 thì nghiệm 휉 = 2 Nếu ( ) ̸= 0 thì chọn [ , ] hoặc [ , ] mà giá trị hàm tại hai đầu trái dấu và ký hiệu là [ 1, 1]. Đối với [ 1, 1] lại tiến hành như [ , ]. Ưu điểm của phương pháp chia đôi là thuật toán đơn giản do đó dễ lập trình. Tuy nhiên do phương pháp chia đôi sử dụng rất ít thông tin về hàm ( ) nên tốc độ hội tụ khá chậm và chỉ sử dụng để giải sơ bộ phương trình.
  43. Matlab trong Giải tích số 24/57 Giải gần đúng phương trình Phương pháp dây cung Phương pháp dây cung Giả sử ( ) liên tục trên đoạn [ , ] và thỏa mãn ( ). ( ) 0. Khi đó, thay vì chia đôi ( ) đoạn [ , ], ta chia theo tỷ lệ − và thu được nghiệm gần đúng ( ) 1 = + ℎ1 trong đó − ( ) ℎ = ( − ). 1 − ( ) + ( ) Tiếp theo dùng cách trên với một trong hai đoạn [ , 1] hay [ 1, ] mà giá trị hàm tại hai đầu trái dấu ta thu được nghiệm gần đúng 2. Công thức lặp của phương pháp dây cung: − 푛−1 푛 = 푛−1 − . ( 푛−1) , ( ) − ( 푛−1) trong đó 0 = (hoặc 0 = ) thì = (hoặc ).
  44. Matlab trong Giải tích số 24/57 Giải gần đúng phương trình Phương pháp dây cung Phương pháp dây cung Giả sử ( ) liên tục trên đoạn [ , ] và thỏa mãn ( ). ( ) 0. Khi đó, thay vì chia đôi ( ) đoạn [ , ], ta chia theo tỷ lệ − và thu được nghiệm gần đúng ( ) 1 = + ℎ1 trong đó − ( ) ℎ = ( − ). 1 − ( ) + ( ) Tiếp theo dùng cách trên với một trong hai đoạn [ , 1] hay [ 1, ] mà giá trị hàm tại hai đầu trái dấu ta thu được nghiệm gần đúng 2. Công thức lặp của phương pháp dây cung: − 푛−1 푛 = 푛−1 − . ( 푛−1) , ( ) − ( 푛−1) trong đó 0 = (hoặc 0 = ) thì = (hoặc ).
  45. Matlab trong Giải tích số 24/57 Giải gần đúng phương trình Phương pháp dây cung Phương pháp dây cung Giả sử ( ) liên tục trên đoạn [ , ] và thỏa mãn ( ). ( ) 0. Khi đó, thay vì chia đôi ( ) đoạn [ , ], ta chia theo tỷ lệ − và thu được nghiệm gần đúng ( ) 1 = + ℎ1 trong đó − ( ) ℎ = ( − ). 1 − ( ) + ( ) Tiếp theo dùng cách trên với một trong hai đoạn [ , 1] hay [ 1, ] mà giá trị hàm tại hai đầu trái dấu ta thu được nghiệm gần đúng 2. Công thức lặp của phương pháp dây cung: − 푛−1 푛 = 푛−1 − . ( 푛−1) , ( ) − ( 푛−1) trong đó 0 = (hoặc 0 = ) thì = (hoặc ).
  46. Matlab trong Giải tích số 25/57 Giải gần đúng phương trình Phương pháp Newton - Raphson Phương pháp Newton - Raphson Phương pháp Newton (còn gọi là phương pháp tiếp tuyến) được dùng nhiều vì ′ nó hội tụ nhanh. Tuy nhiên phương pháp này đòi hỏi phải tính đạo hàm ( ). Công thức Newton - raphson được suy từ khai triển Taylor của ( ) trong lân cận : ′ 2 ( 푖+1) = ( 푖) + ( 푖)( 푖+1 − 푖) + ( 푖+1 − 푖) . Nếu 푖+1 là nghiệm của phương trình ( ) = 0 thì ta có ′ 2 0 = ( 푖) + ( 푖)( 푖+1 − 푖) + ( 푖+1 − 푖) . Giả sử 푖 gần với 푖+1, ta có thể bỏ qua số hạng cuối và thu được công thức Newton - Raphson: ( 푖) 푖+1 = 푖 − ′ (2.2) ( 푖)
  47. Matlab trong Giải tích số 26/57 Giải gần đúng phương trình Phương pháp Newton - Raphson Phương pháp Newton - Raphson Thuật toán được tóm lược như sau 1 Cho 0 ( ) 2 Tính Δ = − ′ ( ) 3 Cho = + Δ 4 Lặp lại bước 2 và bước 3 cho đến khi |Δ | < 휀
  48. Matlab trong Giải tích số 26/57 Giải gần đúng phương trình Phương pháp Newton - Raphson Phương pháp Newton - Raphson Thuật toán được tóm lược như sau 1 Cho 0 ( ) 2 Tính Δ = − ′ ( ) 3 Cho = + Δ 4 Lặp lại bước 2 và bước 3 cho đến khi |Δ | < 휀
  49. Matlab trong Giải tích số 26/57 Giải gần đúng phương trình Phương pháp Newton - Raphson Phương pháp Newton - Raphson Thuật toán được tóm lược như sau 1 Cho 0 ( ) 2 Tính Δ = − ′ ( ) 3 Cho = + Δ 4 Lặp lại bước 2 và bước 3 cho đến khi |Δ | < 휀
  50. Matlab trong Giải tích số 26/57 Giải gần đúng phương trình Phương pháp Newton - Raphson Phương pháp Newton - Raphson Thuật toán được tóm lược như sau 1 Cho 0 ( ) 2 Tính Δ = − ′ ( ) 3 Cho = + Δ 4 Lặp lại bước 2 và bước 3 cho đến khi |Δ | < 휀
  51. Matlab trong Giải tích số 27/57 Giải gần đúng hệ phương trình Nội dung 1 Đa thức nội suy Nội suy Lagrange Nội suy Newton Nội suy bằng đa thức Chebyshev Nội suy bằng đa thức Hermit (đọc thêm) 2 Giải gần đúng phương trình Phương pháp chia đôi Phương pháp dây cung Phương pháp Newton - Raphson 3 Giải gần đúng hệ phương trình Phương pháp lặp đơn Phương pháp Jacobi Phương pháp Gauss-Seidel Phương pháp Newton giải hệ phương trình phi tuyến 4 Giải gần đúng phương trình vi phân thường Phương pháp Euler Phương pháp Euler cải tiến (Modified Euler hay RK-2) Phương pháp Runge-Kutta Hệ phương trình vi phân thường và phương trình vi phân cấp cao
  52. Matlab trong Giải tích số 28/57 Giải gần đúng hệ phương trình Phương pháp lặp đơn Phương pháp lặp đơn Xét hệ phương trình đại số tuyến tính: = , (3.3) 푛×푛 푛 trong đó ∈ R là ma trận cấp 푛 × 푛, ∈ R là vector cho trước, còn 푛 ∈ R là vector nghiệm cần tìm. Để giải lặp hệ (4.1) ta biến đổi nó về dạng thuận tiện cho phép lặp = + , (3.4)
  53. Matlab trong Giải tích số 28/57 Giải gần đúng hệ phương trình Phương pháp lặp đơn Phương pháp lặp đơn Xét hệ phương trình đại số tuyến tính: = , (3.3) 푛×푛 푛 trong đó ∈ R là ma trận cấp 푛 × 푛, ∈ R là vector cho trước, còn 푛 ∈ R là vector nghiệm cần tìm. Để giải lặp hệ (4.1) ta biến đổi nó về dạng thuận tiện cho phép lặp = + , (3.4)
  54. Matlab trong Giải tích số 29/57 Giải gần đúng hệ phương trình Phương pháp lặp đơn Phương pháp lặp đơn Dựa vào nguyên lý ánh xạ co ta có kết quả sau: Định lý 3.1 Giả sử ‖ ‖ < 1. Khi đó mọi dãy lặp ( +1) = ( ) + ( ≥ 0) (3.5) (0) 푛 * trong đó ∈ R bất kỳ cho trước, đều hội tụ đến nghiệm duy nhất của phương trình (3.4). Hơn nữa ta có các đánh giá sai số: ⃦ ⃦ ‖ ‖ ⃦ ⃦ ⃦ ( ) − (0)⃦ ≤ ⃦ (1) − (0)⃦ , ⃦ ⃦ 1 − ‖ ‖ ⃦ ⃦ ⃦ ⃦ ‖ ‖ ⃦ ⃦ ⃦ ( ) − *⃦ ≤ ⃦ ( ) − ( −1)⃦ . ⃦ ⃦ 1 − ‖ ‖ ⃦ ⃦
  55. Matlab trong Giải tích số 30/57 Giải gần đúng hệ phương trình Phương pháp Jacobi Phương pháp Jacobi Định nghĩa 3.1 푛 Ma trận = ( 푖푗 )1 được gọi là ma trận (đường) chéo trội (strictly diagonally dominant) nếu 푛 ∑︁ | 푖푖| > | 푖푗 | với 푖 = 1, 2, . . . , 푛. (3.6) 푗=1 푗̸=푖 Định lý 3.2 Giả sử là ma trận chéo trội. Khi đó hệ (4.1) có thể được biến đổi về dạng phương trình (3.4) với ⎧ 0, 푗 = 푖 푖 푛 ⎨ 푖 = , = ( 푖푗 )1 , 푖푗 = 푖푗 . 푖푖 − , 푗 ̸= 푖 ⎩ 푖푖
  56. Matlab trong Giải tích số 30/57 Giải gần đúng hệ phương trình Phương pháp Jacobi Phương pháp Jacobi Định nghĩa 3.1 푛 Ma trận = ( 푖푗 )1 được gọi là ma trận (đường) chéo trội (strictly diagonally dominant) nếu 푛 ∑︁ | 푖푖| > | 푖푗 | với 푖 = 1, 2, . . . , 푛. (3.6) 푗=1 푗̸=푖 Định lý 3.2 Giả sử là ma trận chéo trội. Khi đó hệ (4.1) có thể được biến đổi về dạng phương trình (3.4) với ⎧ 0, 푗 = 푖 푖 푛 ⎨ 푖 = , = ( 푖푗 )1 , 푖푗 = 푖푗 . 푖푖 − , 푗 ̸= 푖 ⎩ 푖푖
  57. Matlab trong Giải tích số 31/57 Giải gần đúng hệ phương trình Phương pháp Jacobi Phương pháp Jacobi Dễ thấy ∑︀ 푛 | 푖푗 | ∑︁ 푗̸=푖 ‖ ‖∞ = max | 푖푗 | = max < 1. 1≤푖≤푛 1≤푖≤푛 | 푖푖| 푗=1 Công thức lặp Jacobi ( ) (︁ ( ) ( ) ( ))︁ Giả sử phần tử thứ của dãy lặp là = 1 , 2 , . . . , 푛 . Khi đó phần ( +1) (︁ ( +1) ( +1) ( +1))︁ tử kế tiếp = 1 , 2 , . . . , 푛 sẽ được tính theo công thức ( +1) 1 (︁ ( ) ( ) ( ) ( ))︁ 푖 = 푖 − 푖1 1 − · · · − 푖,푖−1 푖−1 − 푖,푖+1 푖+1 − · · · − 푖푛 푛 , 푖푖 푖 = 1, . . . , 푛. (3.7)
  58. Matlab trong Giải tích số 32/57 Giải gần đúng hệ phương trình Phương pháp Gauss-Seidel Phương pháp Gauss-Seidel Công thức lặp Gauss-Seidel ( ) (︁ ( ) ( ) ( ))︁ Giả sử phần tử thứ của dãy lặp là = 1 , 2 , . . . , 푛 . Khi đó phần ( +1) (︁ ( +1) ( +1) ( +1))︁ tử kế tiếp = 1 , 2 , . . . , 푛 sẽ được tính theo công thức ( +1) 1 (︁ ( +1) ( +1) ( ) ( ))︁ 푖 = 푖 − 푖1 1 − · · · − 푖,푖−1 푖−1 − 푖,푖+1 푖+1 − · · · − 푖푛 푛 , 푖푖 푖 = 1, . . . , 푛. (3.8)
  59. Matlab trong Giải tích số 33/57 Giải gần đúng hệ phương trình Phương pháp Newton giải hệ phương trình phi tuyến Phương pháp Newton giải hệ phương trình phi tuyến Để đơn giản, ta chỉ xét trong trường hợp hai chiều. Các kết quả này có thể dễ dàng mở rộng cho các trường hợp nhiều chiều hơn. Mục tiêu Tìm phương pháp giải hệ phương trình phi tuyến 1( , ) = 0 2( , ) = 0, (3.9) trong đó 1, 2 là các hàm (phi tuyến) phụ thuộc vào hai biến , . Định nghĩa 3.2 Ma trận Jacobi của các hàm 1( , ) và 2( , ) được xác định bởi ⎡ 휕 휕 ⎤ 1 ( , ) 1 ( , ) 퐽( , ) = ⎢ 휕 휕 ⎥ (3.10) ⎣ 휕 2 휕 2 ⎦ ( , ) ( , ) 휕 휕
  60. Matlab trong Giải tích số 34/57 Giải gần đúng hệ phương trình Phương pháp Newton giải hệ phương trình phi tuyến Phương pháp Newton giải hệ phương trình phi tuyến Đạo hàm mở rộng Giả sử các hàm = 1( , ) và 푣 = 2( , ) đã biết các giá trị tại điểm ( 0, 0). Ta muốn dự đoán giá trị của chúng tại điểm lân cận ( , ). Khi đó theo các công thức vi phân toàn phần ta có 휕 휕 = 1 ( , ) + 1 ( , ) 휕 0 0 휕 0 0 휕 휕 푣 = 2 ( , ) + 2 ( , ) 휕 0 0 휕 0 0 hoặc có thể viết dưới dạng ma trận [︂ ]︂ [︂ ]︂ 퐹 = = 퐽 ( , ) = 퐽 ( , ) . (3.11) 푣 0 0 0 0
  61. Matlab trong Giải tích số 35/57 Giải gần đúng hệ phương trình Phương pháp Newton giải hệ phương trình phi tuyến Phương pháp Newton giải hệ phương trình phi tuyến Định nghĩa 3.3 Điểm ( , 푞) được gọi là điểm bất động của hai phương trình = 1( , ) (3.12a) = 2( , ) (3.12b) nếu = 1( , 푞) và 푞 = 1( , 푞). Định nghĩa 3.4 Phép lặp điểm bất động của phương trình (3.12a), (3.12b) được xác định bởi +1 = 1 ( , 푞 ) (3.13a) 푞 +1 = 2 ( , 푞 ) , ≥ 0. (3.13b)
  62. Matlab trong Giải tích số 36/57 Giải gần đúng hệ phương trình Phương pháp Newton giải hệ phương trình phi tuyến Phương pháp Newton giải hệ phương trình phi tuyến Định lý 3.3 Giả sử các hàm trong hệ phương trình (3.12a), (3.12b) và các đạo hàm riêng của chúng liên tục trong lân cận của điểm bất động ( , 푞). Khi đó nếu điểm xuất phát ( 0, 푞0) được chọn ”đủ gần” ( , 푞), đồng thời ⃒ ⃒ ⃒ ⃒ ⃒ 휕 1 ⃒ ⃒ 휕 1 ⃒ ⃒ ( , 푞)⃒ + ⃒ ( , 푞)⃒ < 1 ⃒ 휕 ⃒ ⃒ 휕 ⃒ ⃒ ⃒ ⃒ ⃒ ⃒ 휕 2 ⃒ ⃒ 휕 2 ⃒ ⃒ ( , 푞)⃒ + ⃒ ( , 푞)⃒ < 1 ⃒ 휕 ⃒ ⃒ 휕 ⃒ thì dãy lặp (3.13a), (3.13b) hội tụ đến điểm bất động ( , 푞) của hệ (3.12a), (3.12b).
  63. Matlab trong Giải tích số 37/57 Giải gần đúng hệ phương trình Phương pháp Newton giải hệ phương trình phi tuyến Phương pháp Newton giải hệ phương trình phi tuyến Xét các hàm = 1( , ) và 푣 = 2( , ). Giả sử các hàm 1, 2 có các đạo hàm riêng liên tục trong lân cận của điểm ( 0, 0). Khi đó trong lân cận này ta sẽ có các xấp xỉ tuyến tính: 휕 휕 − = 1 ( , )( − ) + 1 ( , )( − ) 0 휕 0 0 0 휕 0 0 0 휕 휕 − = 2 ( , )( − ) + 2 ( , )( − ) , 0 휕 0 0 0 휕 0 0 0 hay có thể viết dưới dạng ⎡ 휕 1 휕 1 ⎤ [︂ ]︂ ( 0, 0) ( 0, 0) [︂ ]︂ − 0 휕 휕 − 0 Δ퐹 = = ⎢ ⎥ = 퐽 ( 0, 0)Δ . ⎣ 휕 2 휕 2 ⎦ 푣 − 푣0 ( , ) ( , ) − 0 휕 0 0 휕 0 0 (3.14)
  64. Matlab trong Giải tích số 37/57 Giải gần đúng hệ phương trình Phương pháp Newton giải hệ phương trình phi tuyến Phương pháp Newton giải hệ phương trình phi tuyến Xét các hàm = 1( , ) và 푣 = 2( , ). Giả sử các hàm 1, 2 có các đạo hàm riêng liên tục trong lân cận của điểm ( 0, 0). Khi đó trong lân cận này ta sẽ có các xấp xỉ tuyến tính: 휕 휕 − = 1 ( , )( − ) + 1 ( , )( − ) 0 휕 0 0 0 휕 0 0 0 휕 휕 − = 2 ( , )( − ) + 2 ( , )( − ) , 0 휕 0 0 0 휕 0 0 0 hay có thể viết dưới dạng ⎡ 휕 1 휕 1 ⎤ [︂ ]︂ ( 0, 0) ( 0, 0) [︂ ]︂ − 0 휕 휕 − 0 Δ퐹 = = ⎢ ⎥ = 퐽 ( 0, 0)Δ . ⎣ 휕 2 휕 2 ⎦ 푣 − 푣0 ( , ) ( , ) − 0 휕 0 0 휕 0 0 (3.14)
  65. Matlab trong Giải tích số 38/57 Giải gần đúng hệ phương trình Phương pháp Newton giải hệ phương trình phi tuyến Phương pháp Newton giải hệ phương trình phi tuyến Ta sẽ dùng (3.14) để xây dựng dãy lặp Newton giải hệ 0 = 1( , ) 0 = 2( , ). (3.15) Giả sử ( , 푞) là nghiệm của (3.15), tức là 0 = 1( , 푞) và 0 = 2( , 푞). Xét sự biến thiên của các hàm , 푣 tại điểm ( 0, 푞0): Δ = − 0 Δ = − 0 Δ푣 = 푣 − 0 Δ푞 = − 푞0 Đặt ( , ) = ( , 푞) trong (3.15) ta thấy ( , 푣) = (0, 0), do đó − 0 = 1( , 푞) − 1 ( 0, 푞0) = − 1( 0, 푞0) 푣 − 푣0 = 2( , 푞) − 2 ( 0, 푞0) = − 2( 0, 푞0)
  66. Matlab trong Giải tích số 39/57 Giải gần đúng hệ phương trình Phương pháp Newton giải hệ phương trình phi tuyến Phương pháp Newton giải hệ phương trình phi tuyến Từ đó ta thu được ⎡ 휕 휕 ⎤ 1 ( , ) 1 ( , ) 0 0 0 0 [︂Δ ]︂ [︂ ( , 푞 )]︂ ⎢ 휕 휕 ⎥ = − 1 0 0 ⎣ 휕 2 휕 2 ⎦ ( , ) ( , ) Δ푞 2 ( 0, 푞0) 휕 0 0 휕 0 0 hay 퐽 ( 0, 푞0)Δ푃 = −퐹 ( 0, 푞0) . (3.16) Nếu ma trận 퐽 ( 0, 푞0) là không suy biến thì ta có −1 Δ푃 = − [퐽 ( 0, 푞0)] 퐹 ( 0, 푞0) . (3.17)
  67. Matlab trong Giải tích số 40/57 Giải gần đúng hệ phương trình Phương pháp Newton giải hệ phương trình phi tuyến Phương pháp Newton giải hệ phương trình phi tuyến Lược đồ lặp Newton Giả sử ở bước thứ ta đã có 푃 . [︂ ]︂ 1 ( , 푞 ) Bước 1: Tính giá trị của vector hàm 퐹 (푃 ) = ; 2 ( , 푞 ) Bước 2: Tính giá trị của ma trận Jacobi ⎡ 휕 휕 ⎤ 1 ( , 푞 ) 1 ( , 푞 ) 휕 휕 퐽 (푃 ) = ⎢ ⎥; ⎣ 휕 2 휕 2 ⎦ ( , 푞 ) ( , 푞 ) 휕 휕 Bước 3: Giải hệ đại số tuyến tính với (ẩn Δ푃 ): 퐽 (푃 )Δ푃 = −퐹 (푃 ); Bước 4: Tính xấp xỉ tiếp theo theo công thức: 푃 +1 = 푃 + Δ푃 ; Lặp lại quá trình trên.
  68. Matlab trong Giải tích số 40/57 Giải gần đúng hệ phương trình Phương pháp Newton giải hệ phương trình phi tuyến Phương pháp Newton giải hệ phương trình phi tuyến Lược đồ lặp Newton Giả sử ở bước thứ ta đã có 푃 . [︂ ]︂ 1 ( , 푞 ) Bước 1: Tính giá trị của vector hàm 퐹 (푃 ) = ; 2 ( , 푞 ) Bước 2: Tính giá trị của ma trận Jacobi ⎡ 휕 휕 ⎤ 1 ( , 푞 ) 1 ( , 푞 ) 휕 휕 퐽 (푃 ) = ⎢ ⎥; ⎣ 휕 2 휕 2 ⎦ ( , 푞 ) ( , 푞 ) 휕 휕 Bước 3: Giải hệ đại số tuyến tính với (ẩn Δ푃 ): 퐽 (푃 )Δ푃 = −퐹 (푃 ); Bước 4: Tính xấp xỉ tiếp theo theo công thức: 푃 +1 = 푃 + Δ푃 ; Lặp lại quá trình trên.
  69. Matlab trong Giải tích số 40/57 Giải gần đúng hệ phương trình Phương pháp Newton giải hệ phương trình phi tuyến Phương pháp Newton giải hệ phương trình phi tuyến Lược đồ lặp Newton Giả sử ở bước thứ ta đã có 푃 . [︂ ]︂ 1 ( , 푞 ) Bước 1: Tính giá trị của vector hàm 퐹 (푃 ) = ; 2 ( , 푞 ) Bước 2: Tính giá trị của ma trận Jacobi ⎡ 휕 휕 ⎤ 1 ( , 푞 ) 1 ( , 푞 ) 휕 휕 퐽 (푃 ) = ⎢ ⎥; ⎣ 휕 2 휕 2 ⎦ ( , 푞 ) ( , 푞 ) 휕 휕 Bước 3: Giải hệ đại số tuyến tính với (ẩn Δ푃 ): 퐽 (푃 )Δ푃 = −퐹 (푃 ); Bước 4: Tính xấp xỉ tiếp theo theo công thức: 푃 +1 = 푃 + Δ푃 ; Lặp lại quá trình trên.
  70. Matlab trong Giải tích số 40/57 Giải gần đúng hệ phương trình Phương pháp Newton giải hệ phương trình phi tuyến Phương pháp Newton giải hệ phương trình phi tuyến Lược đồ lặp Newton Giả sử ở bước thứ ta đã có 푃 . [︂ ]︂ 1 ( , 푞 ) Bước 1: Tính giá trị của vector hàm 퐹 (푃 ) = ; 2 ( , 푞 ) Bước 2: Tính giá trị của ma trận Jacobi ⎡ 휕 휕 ⎤ 1 ( , 푞 ) 1 ( , 푞 ) 휕 휕 퐽 (푃 ) = ⎢ ⎥; ⎣ 휕 2 휕 2 ⎦ ( , 푞 ) ( , 푞 ) 휕 휕 Bước 3: Giải hệ đại số tuyến tính với (ẩn Δ푃 ): 퐽 (푃 )Δ푃 = −퐹 (푃 ); Bước 4: Tính xấp xỉ tiếp theo theo công thức: 푃 +1 = 푃 + Δ푃 ; Lặp lại quá trình trên.
  71. Matlab trong Giải tích số 40/57 Giải gần đúng hệ phương trình Phương pháp Newton giải hệ phương trình phi tuyến Phương pháp Newton giải hệ phương trình phi tuyến Lược đồ lặp Newton Giả sử ở bước thứ ta đã có 푃 . [︂ ]︂ 1 ( , 푞 ) Bước 1: Tính giá trị của vector hàm 퐹 (푃 ) = ; 2 ( , 푞 ) Bước 2: Tính giá trị của ma trận Jacobi ⎡ 휕 휕 ⎤ 1 ( , 푞 ) 1 ( , 푞 ) 휕 휕 퐽 (푃 ) = ⎢ ⎥; ⎣ 휕 2 휕 2 ⎦ ( , 푞 ) ( , 푞 ) 휕 휕 Bước 3: Giải hệ đại số tuyến tính với (ẩn Δ푃 ): 퐽 (푃 )Δ푃 = −퐹 (푃 ); Bước 4: Tính xấp xỉ tiếp theo theo công thức: 푃 +1 = 푃 + Δ푃 ; Lặp lại quá trình trên.
  72. Matlab trong Giải tích số 41/57 Giải gần đúng phương trình vi phân thường Nội dung 1 Đa thức nội suy Nội suy Lagrange Nội suy Newton Nội suy bằng đa thức Chebyshev Nội suy bằng đa thức Hermit (đọc thêm) 2 Giải gần đúng phương trình Phương pháp chia đôi Phương pháp dây cung Phương pháp Newton - Raphson 3 Giải gần đúng hệ phương trình Phương pháp lặp đơn Phương pháp Jacobi Phương pháp Gauss-Seidel Phương pháp Newton giải hệ phương trình phi tuyến 4 Giải gần đúng phương trình vi phân thường Phương pháp Euler Phương pháp Euler cải tiến (Modified Euler hay RK-2) Phương pháp Runge-Kutta Hệ phương trình vi phân thường và phương trình vi phân cấp cao
  73. Matlab trong Giải tích số 42/57 Giải gần đúng phương trình vi phân thường Bài toán giá trị ban đầu Bài toán 4.1 Tìm hàm số = ( ) thỏa mãn điều kiện ˙ = ( , ), 0 ≤ ≤ ; ( 0) = 0. Người ta chứng minh được rằng bài toán trên có duy nhất nghiệm nếu thỏa mãn điều kiện Lipschitz theo đối : | ( , 1) − ( , 2)| ≤ 퐿 | 1 − 2| với 퐿 là một hằng số dương.
  74. Matlab trong Giải tích số 42/57 Giải gần đúng phương trình vi phân thường Bài toán giá trị ban đầu Bài toán 4.1 Tìm hàm số = ( ) thỏa mãn điều kiện ˙ = ( , ), 0 ≤ ≤ ; ( 0) = 0. Người ta chứng minh được rằng bài toán trên có duy nhất nghiệm nếu thỏa mãn điều kiện Lipschitz theo đối : | ( , 1) − ( , 2)| ≤ 퐿 | 1 − 2| với 퐿 là một hằng số dương.
  75. Matlab trong Giải tích số 43/57 Giải gần đúng phương trình vi phân thường Bài toán giá trị ban đầu Mục tiêu Tìm nghiệm bằng số của bài toán (4.1) tại điểm 1 = 0 + ℎ ∈ [ 0, ] với bước ℎ > 0. Một số phương pháp một bước thông dụng Phương pháp Euler Phương pháp Euler cải tiến Phương pháp Runge - Kutta (RK)
  76. Matlab trong Giải tích số 43/57 Giải gần đúng phương trình vi phân thường Bài toán giá trị ban đầu Mục tiêu Tìm nghiệm bằng số của bài toán (4.1) tại điểm 1 = 0 + ℎ ∈ [ 0, ] với bước ℎ > 0. Một số phương pháp một bước thông dụng Phương pháp Euler Phương pháp Euler cải tiến Phương pháp Runge - Kutta (RK)
  77. Matlab trong Giải tích số 43/57 Giải gần đúng phương trình vi phân thường Bài toán giá trị ban đầu Mục tiêu Tìm nghiệm bằng số của bài toán (4.1) tại điểm 1 = 0 + ℎ ∈ [ 0, ] với bước ℎ > 0. Một số phương pháp một bước thông dụng Phương pháp Euler Phương pháp Euler cải tiến Phương pháp Runge - Kutta (RK)
  78. Matlab trong Giải tích số 43/57 Giải gần đúng phương trình vi phân thường Bài toán giá trị ban đầu Mục tiêu Tìm nghiệm bằng số của bài toán (4.1) tại điểm 1 = 0 + ℎ ∈ [ 0, ] với bước ℎ > 0. Một số phương pháp một bước thông dụng Phương pháp Euler Phương pháp Euler cải tiến Phương pháp Runge - Kutta (RK)
  79. Matlab trong Giải tích số 44/57 Giải gần đúng phương trình vi phân thường Phương pháp Euler Phương pháp Euler (RK-1) Công thức Euler ( + ℎ) = ( ) + ℎ ( , ). (4.1) Công thức lặp − Nếu đoạn [ , ] được chia thành đoạn có độ dài ℎ = 0 bởi các điểm 0 chia 0 < 1 < 2 < ··· thì theo công thức Euler ta tính được giá trị gần đúng của nghiệm tại các điểm 푖 = 0 + 푖ℎ, 푖 = 1, 2, theo các công thức: 0 = ( 0); 1 = 0 + ℎ ( 0, 0); ; 푖 = 푖−1 + ℎ ( 푖−1, 푖−1);
  80. Matlab trong Giải tích số 44/57 Giải gần đúng phương trình vi phân thường Phương pháp Euler Phương pháp Euler (RK-1) Công thức Euler ( + ℎ) = ( ) + ℎ ( , ). (4.1) Công thức lặp − Nếu đoạn [ , ] được chia thành đoạn có độ dài ℎ = 0 bởi các điểm 0 chia 0 < 1 < 2 < ··· thì theo công thức Euler ta tính được giá trị gần đúng của nghiệm tại các điểm 푖 = 0 + 푖ℎ, 푖 = 1, 2, theo các công thức: 0 = ( 0); 1 = 0 + ℎ ( 0, 0); ; 푖 = 푖−1 + ℎ ( 푖−1, 푖−1);
  81. Matlab trong Giải tích số 44/57 Giải gần đúng phương trình vi phân thường Phương pháp Euler Phương pháp Euler (RK-1) Công thức Euler ( + ℎ) = ( ) + ℎ ( , ). (4.1) Công thức lặp − Nếu đoạn [ , ] được chia thành đoạn có độ dài ℎ = 0 bởi các điểm 0 chia 0 < 1 < 2 < ··· thì theo công thức Euler ta tính được giá trị gần đúng của nghiệm tại các điểm 푖 = 0 + 푖ℎ, 푖 = 1, 2, theo các công thức: 0 = ( 0); 1 = 0 + ℎ ( 0, 0); ; 푖 = 푖−1 + ℎ ( 푖−1, 푖−1);
  82. Matlab trong Giải tích số 44/57 Giải gần đúng phương trình vi phân thường Phương pháp Euler Phương pháp Euler (RK-1) Công thức Euler ( + ℎ) = ( ) + ℎ ( , ). (4.1) Công thức lặp − Nếu đoạn [ , ] được chia thành đoạn có độ dài ℎ = 0 bởi các điểm 0 chia 0 < 1 < 2 < ··· thì theo công thức Euler ta tính được giá trị gần đúng của nghiệm tại các điểm 푖 = 0 + 푖ℎ, 푖 = 1, 2, theo các công thức: 0 = ( 0); 1 = 0 + ℎ ( 0, 0); ; 푖 = 푖−1 + ℎ ( 푖−1, 푖−1);
  83. Matlab trong Giải tích số 44/57 Giải gần đúng phương trình vi phân thường Phương pháp Euler Phương pháp Euler (RK-1) Công thức Euler ( + ℎ) = ( ) + ℎ ( , ). (4.1) Công thức lặp − Nếu đoạn [ , ] được chia thành đoạn có độ dài ℎ = 0 bởi các điểm 0 chia 0 < 1 < 2 < ··· thì theo công thức Euler ta tính được giá trị gần đúng của nghiệm tại các điểm 푖 = 0 + 푖ℎ, 푖 = 1, 2, theo các công thức: 0 = ( 0); 1 = 0 + ℎ ( 0, 0); ; 푖 = 푖−1 + ℎ ( 푖−1, 푖−1);
  84. Matlab trong Giải tích số 44/57 Giải gần đúng phương trình vi phân thường Phương pháp Euler Phương pháp Euler (RK-1) Công thức Euler ( + ℎ) = ( ) + ℎ ( , ). (4.1) Công thức lặp − Nếu đoạn [ , ] được chia thành đoạn có độ dài ℎ = 0 bởi các điểm 0 chia 0 < 1 < 2 < ··· thì theo công thức Euler ta tính được giá trị gần đúng của nghiệm tại các điểm 푖 = 0 + 푖ℎ, 푖 = 1, 2, theo các công thức: 0 = ( 0); 1 = 0 + ℎ ( 0, 0); ; 푖 = 푖−1 + ℎ ( 푖−1, 푖−1);
  85. Matlab trong Giải tích số 45/57 Giải gần đúng phương trình vi phân thường Phương pháp Euler cải tiến (Modified Euler hay RK-2) Phương pháp Euler cải tiến Công thức Euler cải tiến = ( ) + ℎ ( , ) ℎ ( + ℎ) = [ ( , ) + ( + ℎ, )] . (4.2) 2 Công thức lặp Với 푖 = 1, 2, : = 푖−1 + ℎ ( 푖−1, 푖−1); 푖 = 푖−1 + ℎ [ ( 푖−1, 푖−1) + ( 푖−1 + ℎ, )] (4.3)
  86. Matlab trong Giải tích số 45/57 Giải gần đúng phương trình vi phân thường Phương pháp Euler cải tiến (Modified Euler hay RK-2) Phương pháp Euler cải tiến Công thức Euler cải tiến = ( ) + ℎ ( , ) ℎ ( + ℎ) = [ ( , ) + ( + ℎ, )] . (4.2) 2 Công thức lặp Với 푖 = 1, 2, : = 푖−1 + ℎ ( 푖−1, 푖−1); 푖 = 푖−1 + ℎ [ ( 푖−1, 푖−1) + ( 푖−1 + ℎ, )] (4.3)
  87. 2 = ( 푛 + 2ℎ, 푛 + ℎ 21 1) 3 = ( 푛 + 3ℎ, 푛 + ℎ ( 31 1 + 32 2)) . . (︃ 푠−1 )︃ ∑︁ 푠 = 푛 + 푠ℎ, 푛 + ℎ 푠푖 푖 푖=1 Matlab trong Giải tích số 46/57 Giải gần đúng phương trình vi phân thường Phương pháp Runge-Kutta Phương pháp Runge-Kutta Họ các phương pháp Runge-Kutta (RK) hiển s-nấc được xác định bởi công thức 푠 ∑︁ 푛+1 = 푛 + ℎ 푖 푖, (4.4) 푖=1 trong đó 1 = ( 푛, 푛)
  88. 3 = ( 푛 + 3ℎ, 푛 + ℎ ( 31 1 + 32 2)) . . (︃ 푠−1 )︃ ∑︁ 푠 = 푛 + 푠ℎ, 푛 + ℎ 푠푖 푖 푖=1 Matlab trong Giải tích số 46/57 Giải gần đúng phương trình vi phân thường Phương pháp Runge-Kutta Phương pháp Runge-Kutta Họ các phương pháp Runge-Kutta (RK) hiển s-nấc được xác định bởi công thức 푠 ∑︁ 푛+1 = 푛 + ℎ 푖 푖, (4.4) 푖=1 trong đó 1 = ( 푛, 푛) 2 = ( 푛 + 2ℎ, 푛 + ℎ 21 1)
  89. . . (︃ 푠−1 )︃ ∑︁ 푠 = 푛 + 푠ℎ, 푛 + ℎ 푠푖 푖 푖=1 Matlab trong Giải tích số 46/57 Giải gần đúng phương trình vi phân thường Phương pháp Runge-Kutta Phương pháp Runge-Kutta Họ các phương pháp Runge-Kutta (RK) hiển s-nấc được xác định bởi công thức 푠 ∑︁ 푛+1 = 푛 + ℎ 푖 푖, (4.4) 푖=1 trong đó 1 = ( 푛, 푛) 2 = ( 푛 + 2ℎ, 푛 + ℎ 21 1) 3 = ( 푛 + 3ℎ, 푛 + ℎ ( 31 1 + 32 2))
  90. (︃ 푠−1 )︃ ∑︁ 푠 = 푛 + 푠ℎ, 푛 + ℎ 푠푖 푖 푖=1 Matlab trong Giải tích số 46/57 Giải gần đúng phương trình vi phân thường Phương pháp Runge-Kutta Phương pháp Runge-Kutta Họ các phương pháp Runge-Kutta (RK) hiển s-nấc được xác định bởi công thức 푠 ∑︁ 푛+1 = 푛 + ℎ 푖 푖, (4.4) 푖=1 trong đó 1 = ( 푛, 푛) 2 = ( 푛 + 2ℎ, 푛 + ℎ 21 1) 3 = ( 푛 + 3ℎ, 푛 + ℎ ( 31 1 + 32 2)) . .
  91. Matlab trong Giải tích số 46/57 Giải gần đúng phương trình vi phân thường Phương pháp Runge-Kutta Phương pháp Runge-Kutta Họ các phương pháp Runge-Kutta (RK) hiển s-nấc được xác định bởi công thức 푠 ∑︁ 푛+1 = 푛 + ℎ 푖 푖, (4.4) 푖=1 trong đó 1 = ( 푛, 푛) 2 = ( 푛 + 2ℎ, 푛 + ℎ 21 1) 3 = ( 푛 + 3ℎ, 푛 + ℎ ( 31 1 + 32 2)) . . (︃ 푠−1 )︃ ∑︁ 푠 = 푛 + 푠ℎ, 푛 + ℎ 푠푖 푖 푖=1
  92. hoặc ở dạng gọn hơn: với = ( 푖푗 ) và 푖푗 = 0 với 푗 ≥ 푖. Matlab trong Giải tích số 47/57 Giải gần đúng phương trình vi phân thường Phương pháp Runge-Kutta Phương pháp Runge-Kutta Bảng Butcher 0 2 21 3 31 32 . . . . . . . . . . . . . 푠 푠1 푠2 ······ 푠,푠−1 1 2 ······ 푠−1 푠
  93. Matlab trong Giải tích số 47/57 Giải gần đúng phương trình vi phân thường Phương pháp Runge-Kutta Phương pháp Runge-Kutta Bảng Butcher 0 hoặc ở dạng gọn hơn: 2 21 3 31 32 với = ( 푖푗 ) và 푖푗 = 0 với 푗 ≥ 푖. . . . . . . . . . . . . . 푠 푠1 푠2 ······ 푠,푠−1 1 2 ······ 푠−1 푠
  94. Matlab trong Giải tích số 48/57 Giải gần đúng phương trình vi phân thường Phương pháp Runge-Kutta Phương pháp Runge-Kutta Các ví dụ Ví dụ 1 Xét trường hợp 푠 = 1. Khi đó 1 = ( 푛, 푛) 푛+1 = 푛 + ℎ 1 ( 푛, 푛) . Mặt khác, áp dụng công thức khai triển Taylor: 2 푛+1 = 푛 + ℎ ˙| 푛 + ··· = 푛 + ℎ ( 푛, 푛) + 풪(ℎ ) =⇒ 1 = 1. Do đó, phương pháp Runge - Kutta một nấc (RK 1) tương đương với phương pháp Euler hiển.
  95. Matlab trong Giải tích số 49/57 Giải gần đúng phương trình vi phân thường Phương pháp Runge-Kutta Phương pháp Runge-Kutta Các ví dụ Ví dụ 2 Với 푠 = 2, phương trình (4.4) tương đương với hệ 1 = ( 푛, 푛) 2 = ( 푛 + 2ℎ, 푛 + ℎ 21 1) , 푛+1 = 푛 + ℎ ( 1 1 + 2 2) . (4.5) Áp dụng khai triển Taylor trong lân cận 푛 ta có ℎ2 2 = + ℎ | + | + 풪 (︀ℎ3)︀ . 푛+1 푛 푛 2 2 푛 Mặt khác, ta biết rằng ˙ = ( , ), do đó 2 ( , ) 휕 ( , ) 휕 ( , ) := = + ( , ) . 2 휕 휕
  96. Matlab trong Giải tích số 50/57 Giải gần đúng phương trình vi phân thường Phương pháp Runge-Kutta Phương pháp Runge-Kutta Các ví dụ Ví dụ 2 (tiếp) Do đó, công thức khai triển Taylor có thể được viết lại như sau ℎ2 (︂ 휕 휕 )︂ − = ℎ ( , ) + + | + 풪 (︀ℎ3)︀ . (4.6) 푛+1 푛 푛 푛 2 휕 휕 ( 푛, 푛) (︀ 3)︀ Mặt khác, số hạng 2 trong công thức RK trên có thể khai triển tới 풪 ℎ bởi 2 = ( 푛 + 2ℎ, 푛 + ℎ 21 1) 휕 휕 = ℎ ( , ) + ℎ | + +ℎ | + 풪 (︀ℎ3)︀ . 푛 푛 2 휕 ( 푛, 푛) 21 휕 ( 푛, 푛) Thay vào phương trình cuối của hệ (4.5) ta thu được 휕 휕 − = ℎ ( + ) ( , )+ℎ2 | +ℎ2 | +풪 (︀ℎ3)︀ . 푛+1 푛 1 2 푛 푛 2 2 휕 ( 푛, 푛) 2 21 휕 ( 푛, 푛) Thay vào (4.6) ta thu được hệ
  97. Matlab trong Giải tích số 51/57 Giải gần đúng phương trình vi phân thường Phương pháp Runge-Kutta Phương pháp Runge-Kutta Các ví dụ Ví dụ 2 (tiếp) 1 + 2 = 1, 1 = , 2 2 2 1 = . 2 21 2 1 1 Cho = 1. Khi đó = , = , = 1. Bảng Butcher tương ứng 2 2 2 1 2 21 0 1 1 . Do đó, công thức Runge-Kutta 2-nấc trong trường hợp này có 1/2 1/2 dạng công thức Heun: ℎ = + ( ( , ) + ( + ℎ, + ℎ ( , ))) , 푛+1 푛 2 푛 푛 푛 푛 푛 푛
  98. Matlab trong Giải tích số 52/57 Giải gần đúng phương trình vi phân thường Phương pháp Runge-Kutta Phương pháp RK4 Công thức RK4 1 = ℎ ( , ); (︂ ℎ )︂ = ℎ + , + 1 ; 2 2 2 (︂ ℎ )︂ = ℎ + , + 2 ; 3 2 2 4 = ℎ ( + ℎ, + 3); 1 = ( ) + ( + 2 + 2 + ) . (4.7) 6 1 2 3 4 0 1/2 1/2 Bảng Butcher tương ứng 1/2 0 1/2 1 0 0 1 1/6 1/3 1/3 1/6
  99. Matlab trong Giải tích số 53/57 Giải gần đúng phương trình vi phân thường Phương pháp Runge-Kutta Phương pháp RK4 Công thức lặp Với 푖 = 1, 2 : 1 = ℎ ( 푖−1, 푖−1); (︂ ℎ )︂ = ℎ + , + 1 ; 2 푖−1 2 푖−1 2 (︂ ℎ )︂ = ℎ + , + 2 ; 3 푖−1 2 푖−1 2 4 = ℎ ( 푖−1 + ℎ, 푖−1 + 3); 1 = + ( + 2 + 2 + ) . (4.8) 푖 푖−1 6 1 2 3 4
  100. Matlab trong Giải tích số 54/57 Giải gần đúng phương trình vi phân thường Hệ phương trình vi phân thường và phương trình vi phân cấp cao Bài toán Bài toán 4.2 Tìm = ( ) và = ( ) là nghiệm của bài toán Cauchy ˙ = ( , , ) ( 0) = 훼 ˙ = ( , , ) ( 0) = 훽 ∈ [ 0, ] (4.9) hoặc Bài toán 4.3 Tìm = ( ) là nghiệm của bài toán Cauchy cấp hai: ¨ = ( , , ) ∈ [ 0, ] ( 0) = 훼; ˙( 0) = 훽. (4.10)
  101. Matlab trong Giải tích số 54/57 Giải gần đúng phương trình vi phân thường Hệ phương trình vi phân thường và phương trình vi phân cấp cao Bài toán Bài toán 4.2 Tìm = ( ) và = ( ) là nghiệm của bài toán Cauchy ˙ = ( , , ) ( 0) = 훼 ˙ = ( , , ) ( 0) = 훽 ∈ [ 0, ] (4.9) hoặc Bài toán 4.3 Tìm = ( ) là nghiệm của bài toán Cauchy cấp hai: ¨ = ( , , ) ∈ [ 0, ] ( 0) = 훼; ˙( 0) = 훽. (4.10)
  102. Matlab trong Giải tích số 55/57 Giải gần đúng phương trình vi phân thường Hệ phương trình vi phân thường và phương trình vi phân cấp cao Bài toán Ta thấy bài toán (4.3) có thể được đưa về bài toán (4.2) bằng cách đặt ˙ = ( ) ( 0) = 훼 ˙ = ( , , ) ( 0) = 훽. (4.11) Vậy ta chỉ xét phương pháp tìm nghiệm của hệ hai phương trình vi phân cấp 1 (4.9).
  103. Matlab trong Giải tích số 56/57 Giải gần đúng phương trình vi phân thường Hệ phương trình vi phân thường và phương trình vi phân cấp cao Phương pháp Runge - Kutta giải hệ (4.9) Chia đoạn [ 0, ] thành 푛 đoạn con bởi các điểm chia 0 < 1 < ··· < 푛. Độ dài mỗi đoạn ℎ = 푖+1 − 푖. Giả sử giá trị của nghiệm tại 푖−1 là 푖−1 và 푖−1 đã biết, ta tìm 푖 ≈ ( 푖) và 푖 ≈ ( 푖) theo các công thức của phương pháp RK4: 1 = ℎ ( 푖−1, 푖−1, 푖−1); 푙1 = ℎ ( 푖−1, 푖−1, 푖−1) (︂ ℎ 푙 )︂ = ℎ + , + 1 , + 1 ; 2 푖−1 2 푖−1 2 푖−1 2 (︂ ℎ 푙 )︂ 푙 = ℎ + , + 1 , + 1 2 푖−1 2 푖−1 2 푖−1 2 (︂ ℎ 푙 )︂ = ℎ + , + 2 , + 2 ; 3 푖−1 2 푖−1 2 푖−1 2 (︂ ℎ 푙 )︂ 푙 = ℎ + , + 2 , + 2 3 푖−1 2 푖−1 2 푖−1 2 4 = ℎ ( 푖−1 + ℎ, 푖−1 + 3, 푖−1 + 푙3); 푙4 = ℎ ( 푖−1 + ℎ, 푖−1 + 3, 푖−1 + 푙3) 1 1 = + ( + 2 + 2 + ); = + (푙 + 2푙 + 2푙 + 푙 )(4.12). 푖 푖−1 6 1 2 3 4 푖 푖−1 6 1 2 3 4
  104. Matlab trong Giải tích số 57/57 Giải gần đúng phương trình vi phân thường Hệ phương trình vi phân thường và phương trình vi phân cấp cao Bài tập thực hành Tìm nghiệm của bài toán Cauchy cấp hai: ¨ + ˙ + = 0 0 ≤ ≤ 0.5 (0) = 0 ˙(0) = 1. (4.13) Bài toán đã cho tương đương với hệ hai phương trình vi phân cấp 1 ˙ = , (0) = 0 ˙ = − − , (0) = 1. (4.14)