Tài liệu Tính ổn định và tính co của các phương pháp runge kutta

ĐẠI HỌC QUỐC GIA HÀ NỘI TRƯỜNG ĐẠI HỌC KHOA HỌC TỰ NHIÊN ———————————— NGUYỄN THỊ HIÊN TÍNH ỔN ĐỊNH VÀ TÍNH CO CỦA CÁC PHƯƠNG PHÁP RUNGE-KUTTA Chuyên ngành: Toán ứng dụng Mã số: 60 46 01 12 LUẬN VĂN THẠC SĨ KHOA HỌC NGƯỜI HƯỚNG DẪN KHOA HỌC: PGS.TS. VŨ HOÀNG LINH Hà Nội-2014 Lời cảm ơn Trước khi trình bày nội dung chính của luận văn, tôi xin cảm ơn Ban chủ nhiệm khoa Toán - Cơ - Tin học cùng toàn thể các thầy giáo, cô giáo trong khoa Toán - Cơ - Tin học, phòng Sau đại học, trường Đại học Khoa học Tự nhiên - Đại học Quốc gia Hà Nội, đã giảng dạy tận tình và tạo điều kiện thuận lợi để tôi hoàn thành tốt luận văn. Đặc biệt, tôi xin gửi lời cảm ơn chân thành nhất tới thầy giáo PGS.TS Vũ Hoàng Linh, người đã trực tiếp hướng dẫn và tận tình chỉ bảo tôi trong suốt quá trình tôi học tập và thực hiện luận văn. Nhân dịp này, tôi cũng xin cảm ơn gia đình đã luôn ủng hộ và động viên trong suốt thời gian tôi học tập. Cuối cùng, tôi xin cảm ơn tất cả các bạn, các anh, các chị trong lớp cao học Toán khóa 2011 - 2013, đặc biệt là các anh chị chuyên ngành Toán ứng dụng khóa 2010 - 2012 và khóa 2011 - 2013 đã tận tình giúp đỡ và động viên tôi trong quá trình học tập. Xin chân thành cảm ơn! Hà Nội, ngày 16 tháng 01 năm 2014 Học viên Nguyễn Thị Hiên Mục lục Mở đầu 3 Bảng ký hiệu 4 1 Các khái niệm cơ bản 5 1.1 Các phương pháp Runge-Kutta . . . . . . . . . . . . . . . . 5 1.2 Xây dựng các phương pháp Runge-Kutta ẩn . . . . . . . . 11 1.3 Áp dụng các phương pháp Runge-Kutta giải bài toán cương 18 1.4 Các loại chuẩn . . . . . . . . . . . . . . . . . . . . . . . . 21 2 Tính co cho bài toán tuyến tính 2.1 Chuẩn Euclid (Định lý von Neumann) . . . . . 2.2 Hàm tăng trưởng sai số với bài toán tuyến tính 2.3 Bài toán với nhiễu phi tuyến nhỏ . . . . . . . . 2.4 Tính co trong k.k∞ và k.k1 . . . . . . . . . . . 2.5 Hệ số ngưỡng . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 Tính ổn định B và tính co 3.1 Điều kiện Lipschitz một phía . . . . . . . . . . . . . . 3.2 Ổn định B và ổn định đại số . . . . . . . . . . . . . . 3.3 Một vài phương pháp Runge-Kutta ẩn ổn định đại số . 3.4 Ổn định AN . . . . . . . . . . . . . . . . . . . . . . . 3.5 Các phương pháp Runge-Kutta khả quy . . . . . . . . 3.6 Định lý về sự tương đương giữa ổn định B và ổn định số với các phương pháp S-bất khả quy . . . . . . . . . 3.7 Hàm tăng trưởng sai số . . . . . . . . . . . . . . . . . 3.8 Tính toán ϕB (x) . . . . . . . . . . . . . . . . . . . . Kết luận . . . . . . . . . . . . . . . . . . . . . . . . . . . Tài liệu tham khảo . . . . . . . . . . . . . . . . . . . . . 2 . . . . . . . . . . . . . . . . . . . . đại . . . . . . . . . . . . . . . 26 29 30 33 37 39 . . . . . 42 42 43 46 48 51 . . . . . 53 56 58 63 64 Mở đầu Trong khoa học và kĩ thuật ta thường gặp rất nhiều bài toán liên quan tới việc giải phương trình vi phân. Có rất nhiều trường hợp nghiệm giải tích của các bài toán này là không thể tìm được. Chính vì vậy các nhà toán học đã tìm kiếm nhiều phương pháp số khác nhau để giải các bài toán trên. Trong các phương pháp số, phương pháp Runge-Kutta có nhiều tính chất ưu việt và được sử dụng rộng rãi. Luận văn trình bày về tính ổn định và tính co của các phương pháp Runge-Kutta. Xuất phát từ điều kiện ổn định tuyệt đối |yn | ≤ |yn−1 | của bài toán y 0 = λy , ta mở rộng đến khái niệm "tính co" khi xét bài toán tuyến tính y 0 = Ay , tiếp đến là các khái niệm tính ổn định B và ổn định đại số khi xét bài toán phi tuyến. Trên cơ sở đó ta có thể lựa chọn ra phương pháp hữu hiệu và phù hợp nhất để giải các bài toán nảy sinh trong thực tế. Nội dung luận văn được tham khảo chính từ tài liệu [2] và [3]. Bố cục của luận văn bao gồm 3 chương: • Chương 1: Các khái niệm cơ bản Luận văn trình bày các khái niệm cơ bản về phương pháp Runge-Kutta, cách xây dựng phương pháp Runge-Kutta ẩn, cùng với các kiến thức bổ trợ cho Chương 2 và Chương 3. • Chương 2: Tính co của bài toán tuyến tính Luận văn trình bày các khái niệm và định lý liên quan đến tính co khi xét bài toán tuyến tính. • Chương 3: Tính ổn định B và tính co Luận văn trình bày khái niệm ổn định B, ổn định đại số, ổn định AN và mối quan hệ giữa các khái niệm ổn định của các phương pháp Runge-Kutta khi xét bài toán phi tuyến. Do thời gian thực hiện luận văn không nhiều nên trong luận văn không tránh khỏi những hạn chế và sai sót. Tác giả mong nhận được sự góp ý và những ý kiến phản biện của quý thầy cô và bạn đọc. 3 Bảng ký hiệu A⊗I Tích tensor. B (p), C (η), D (ζ) Bộ điều kiện về cấp chính xác. C Tập số phức. Cn Không gian vectơ phức n chiều. I Ma trận đơn vị. K (Z) Hàm ổn định với bài toán y 0 = λ (x) y. Pk (x) Đa thức trực giao Legendre. Pkj (z) Xấp xỉ Padé. R (z) Hàm ổn định của phương pháp. R Tập số thực. Rn Không gian vectơ thực n chiều. S Miền ổn định. µ (A) Chuẩn logarit của ma trận A. ν Hằng số Lipschitz một phía. ϕB (x) Hàm tăng trưởng sai số khi xét bài toán phi tuyến. ϕR (x) Hàm tăng trưởng sai số khi xét bài toán tuyến tính. % Hệ số ngưỡng.  bT = (b1 , ..., bs )   Chuyển vị của vectơ b =    b1 .. .    .   bs T 1 = (1, ... , 1) Vectơ cột với các tất cả các thành phần đều bằng 1. 4 Chương 1 Các khái niệm cơ bản Chương này trình bày các khái niệm cơ bản về các phương pháp RungeKutta, sự tồn tại lời giải số của phương pháp, cách xây dựng các phương pháp Runge-Kutta ẩn cùng với các kiến thức bổ trợ cho Chương 2 và Chương 3. Nội dung của chương này chỉ phát biểu các khái niệm và các kết quả phục vụ cho các chương sau. Chứng minh chi tiết của các kết quả trong chương này có thể tham khảo tại [2], [3] và [5]. 1.1 Các phương pháp Runge-Kutta Phương pháp Runge-Kutta tổng quát Phương pháp Runge-Kutta thuộc lớp các phương pháp số một bước, được đưa ra bởi hai nhà toán học người Đức là Carl Runge (1856 - 1927) và Wilhelm Kutta (1867 - 1944). Trước hết ta xét bài toán Cauchy của phương trình vi phân cấp một có dạng y 0 = f (t, y) , y ∈ Rn , f : R × Rn → Rn , y (t0 ) = y0 . (1.1) Định nghĩa 1.1 (xem [5]). Phương pháp Runge-Kutta s nấc cho hệ phương trình vi phân (1.1) có thể viết dưới dạng: Yi = yn−1 + h yn = yn−1 + h s P j=1 s P aij f (tn−1 + cj h, Yj ) i = 1, ..., s (1.2) bi f (tn−1 + ci h, Yi ) i=1 Trong đó Y1 , ..., Ys là các giá trị nấc xấp xỉ y tại ti = tn−1 + ci h (ti là các điểm s P nấc). Bộ hệ số: {ci }si=1 ; {aij }si,j=1 ; {bi }si=1 thỏa mãn bi = 1. Thông thường, ta chọn để có ci = s P i=1 aij (i = 1, ..., s). j=1 5 • Nếu aij = 0 với i ≤ j thì phương pháp là phương pháp Runge-Kutta hiển (ERK). • Nếu aij = 0 với i < j và có ít nhất một aii 6= 0 thì phương pháp là phương pháp Runge-Kutta ẩn đường chéo (DIRK). • Nếu aij = 0 với i < j và aii = γ với i = 1, ..., s thì phương pháp là phương pháp ẩn đường chéo đơn (SDIRK). • Các trường hợp còn lại gọi là phương pháp Runge-Kutta ẩn (IRK). Để dễ dàng hình dung về phương pháp Runge-Kutta, Butcher đã đưa bộ hệ số của phương pháp vào bảng sau: c1 c2 .. . a11 a21 .. . a12 a22 .. . ··· ··· .. . a1s a2s .. . cs as1 b1 as2 b2 ··· ··· ass bs Bảng 1.1: Bảng Butcher. Ví dụ 1.1. Một số công thức ERK cơ bản (a) Euler hiển 0 0 1 (b) Hình thang hiển 0 1 0 1 1 2 (c) Trung điểm hiển 0 1 2 0 0 1 2 0 1 2 0 0 0 1 Bảng 1.2: Một số công thức ERK. Ví dụ 1.2. Một số công thức IRK cơ bản (a) Euler ẩn 1 1 1 (b) Hình thang ẩn 0 1 0 1 2 1 2 0 1 2 1 2 Bảng 1.3: Một số công thức IRK. 6 (c) Trung điểm ẩn 1 2 1 2 1 Sự tồn tại lời giải số của phương pháp Xét công thức (1.2) trong trường hợp n = 1, nếu ta đặt ki = f (t0 + ci h, Yi ) với i = 1, 2, ..., s ta thu được s P ki = f (t0 + ci h, y0 + h aij kj ) j=1 y1 = y0 + s P (1.3) bi k i . i=1 Để xác định lời giải số y1 của phương pháp, trước hết ta cần xác định các giá trị ki từ hệ phương trình chứa ki cho bởi (1.3). Nói chung, đây là hệ phương trình phi tuyến nên trong nhiều trường hợp có thể các ki tồn tại không duy nhất. Do đó, không tồn tại lời giải số của phương pháp. Định lý sau đây cho ta điều kiện để tồn tại lời giải số của phương pháp Runge-Kutta ẩn (1.3). Định lý 1.1 (xem [2]). Cho hàm f : R × Rn → Rn là hàm liên tục và thỏa mãn điều kiện Lipschitz theo biến y với hằng số L. Nếu h< 1 P , L max ti |aij | j thì lời giải số của phương pháp (1.3) tồn tại duy nhất với các giá trị ki được xác định duy nhất từ hệ phương trình cho bởi (1.3), các giá trị này có thể thu được bằng phương pháp lặp Newton. Hơn nữa, nếu f (t, y) là hàm khả vi, liên tục tới cấp p thì các ki (là các hàm theo biến h) cũng khả vi, liên tục cho tới cấp p. Ổn định tuyệt đối và ổn định A Bằng các phương pháp khác nhau ta có thể tìm ra nghiệm số của phương trình vi phân. Tuy nhiên, nghiệm số ta tìm được liệu có tốt không, làm thế nào để có thể đánh giá được nghiệm đó. Để giải quyết vấn đề này, ta cần chỉ ra nghiệm số đó phải có tính chất tốt cho ít nhất một lớp các bài toán. Xét phương trình vi phân thử y 0 = λy, λ là hằng số, λ ∈ C, y(t0 ) = y0 . (1.4) Ở các phần sau, khi xét phương trình thử (1.4) ta luôn giả sử Re (λ) ≤ 0. Trong trường hợp Re (λ) ≤ 0 ta được điều kiện ổn định tuyệt đối là |yn | ≤ |yn−1 | , n = 1, 2, .... (1.5) Giả sử y(t), ye(t) là hai lời giải của (1.1). Từ đó ta có các định nghĩa về sự ổn định và ổn định tiệm cận với nghiệm của phương trình vi phân( 1.1). 7 Định nghĩa 1.2. Nghiệm y (t) của phương trình vi phân( 1.1) gọi là ổn định nếu với mọi ε > 0, ∃δ > 0 sao cho |y (to ) − ye (to )| ≤ δ thì |y (t) − ye (t)| < ε với mọi t ≥ t0 . Định nghĩa 1.3. Nghiệm y (t) của phương trình vi phân( 1.1) gọi là ổn định tiệm cận nếu nó ổn định và thỏa mãn điều kiện lim |y (t) − ye (t)| = 0. t→+∞ Nhận xét 1.1. Đối với hệ tuyến tính, nếu nghiệm tầm thường là ổn định (ổn định tiệm cận) thì mọi nghiệm cũng ổn định (ổn định tiệm cận). Trong trường hợp này, chúng ta nói hệ ổn định (ổn định tiệm cận). Xét bài toán tuyến tính y 0 = Ay , A ∈ Cm×m . (1.6) Định lý 1.2 (xem [5]). Nếu mọi giá trị riêng λ của ma trận A đều thỏa mãn Re (λ) ≤ 0 và các giá trị riêng có phần thực bằng 0 đều là các giá trị riêng đơn thì hệ y 0 = Ay là ổn định. Định lý 1.3 (Điều kiện cần và đủ (xem [5])). Nếu Re (λ) < 0 với mọi giá trị riêng λ của A thì hệ ổn định tiệm cận. Để đi đến khái niệm về hàm ổn định của phương pháp Runge-Kutta, ta áp dụng phương pháp Runge-Kutta với công thức (1.2) cho phương trình thử (1.4) được lời giải số yn = R (z) yn−1 với z = λh. (1.7) Khi đó, để điều kiện ổn định tuyệt đối (1.5) được thỏa mãn thì |R (z)| ≤ 1. (1.8) Định nghĩa 1.4. Hàm R (z) xác định ở (1.7) được gọi là hàm ổn định của phương pháp Runge-Kutta. Tập S = {z ∈ C : |R (z)| ≤ 1} được gọi là miền ổn định tuyệt đối của phương pháp Runge-Kutta. Ví dụ 1.3. • Phương pháp Euler hiển có hàm ổn định R (z) = 1 + z . • Phương pháp Euler hiển có miền ổn định tuyệt đối là hình tròn có bán kính bằng 1, tâm −1. 8 Mệnh đề 1.1 (xem [3]). Phương pháp Runge-Kutta ẩn s nấc với s X gi = y 0 + h aij f (t0 + cj h, gj ), j=1 s X y1 = y0 + h i = 1, 2, ..., s, (1.9a) (1.9b) bj f (t0 + cj h, gj ), j=1 khi áp dụng cho phương trình thử y 0 = λy thì có y1 = R (hλ) y0 với R (z) = 1 + zbT (I − zA)−11 , trong đó bT = (b1 , ..., bs ) , A = (aij )si,j=1 , (1.10) 1 = (1, 1, ..., 1)T . Mệnh đề 1.2 (xem [3]). Hàm ổn định của (1.9) thỏa mãn
T R (z) = STT Phương pháp 1 Euler ẩn 2 Phương pháp θ 3 4 Trung điểm ẩn Hình thang ẩn det I − zA + z11b det (I − zA)  SDIRK cấp 3 (1.11) . R (z) 1 1−z 1 + z (1 − θ) 1 − zθ 1 + z/2 1 − z/2
1 + z (1 − 2γ) + z 2 1/2 − 2γ + γ 2 2 (1 − γz) Bảng 1.4: Hàm ổn định của một số phương pháp Runge-Kutta ẩn. Miền ổn định tuyệt đối của các phương pháp cho trong Bảng 1.4 có thể xem trong [3]. Từ kết quả ở trên ta thấy phương pháp Runge-Kutta ẩn có hàm ổn định R (z) là hàm hữu tỉ với tử số và mẫu số có bậc ≤ s. R (z) = P (z) , Q (z) deg P = k, deg Q = j. (1.12) Định nghĩa 1.5. Một phương pháp mà miền ổn định tuyệt đối thỏa mãn S ⊃ C− = {z ∈ C : Re (z) ≤ 0} , thì phương pháp đó gọi là phương pháp ổn định A (hay A-ổn định). Phương pháp Runge-Kutta với hàm ổn định như ở (1.12) là ổn định A khi và chỉ khi |R (iy)| ≤ 1 với mọi y ∈ R (1.13) và R (z) là giải tích với Re (z) < 0. 9 (1.14) Xấp xỉ Padé của ez Cho hàm f (z) giải tích trong lân cận của điểm 0, hai số nguyên không âm k và j . Khi đó, ta có thể xấp xỉ hàm f (z) bởi hàm hữu tỉ f (z) ≈ Pk (z) . Qj (z) (1.15)
P là đa thức bậc k , Q là đa thức bậc j và sai số của xấp xỉ tới O z k+j+1 . Trong trường hợp j = 0 khi đó xấp xỉ chính là khai triển Taylor của hàm f (z) tại điểm 0. Khi k = 0 thì Q (z)/P (z) là khai triển Taylor của hàm 1/f (z). Tuy nhiên, với hàm bất kỳ và k, j bất kỳ thì không phải lúc nào cũng tồn tại xấp xỉ này. Khi xấp xỉ f (z) =
Pkj (z) + O z k+j+1 Qkj (z) tồn tại, thì bộ số (k , j ) được gọi là xấp xỉ Padé của hàm f . Xấp xỉ Padé của hàm mũ là trường hợp đặc biệt chúng ta cần nghiên cứu, một vài xấp xỉ Padé xấp xỉ các hàm hữu tỉ của các phương pháp quan trọng như phương pháp Gauss, phương pháp Radau, phương pháp Lobatto... Ta sẽ chỉ ra luôn tồn tại xấp xỉ Padé cho các hàm ổn định của các phương pháp này với k , j bất kỳ. Định lý 1.4 (xem [3]). Bộ số (k , j ) là xấp xỉ Padé của ez được cho bởi Rkj (z) = Pkj (z) , Qkj (z) (1.16) trong đó k k (k − 1) z2 k (k − 1) ...1 z k .z + . + ... + . j+k (j + k) (j + k − 1) 2! (j + k) ... (j + 1) k! j (j − 1) z2 (−1)j j (j − 1) ...1 z j j =1− .z + . + ... + . = Pjk (−z) . k+j (k + j) (k + j − 1) 2! (k + j) ... (k + 1) j! Pkj = 1 + Qkj Đặt Ckj = (−1)j k!j! , (k + j)! (k + j + 1)! khi đó Pkj − ez Qkj + Ckj z k+j+1 +
j+1 Ckj z k+j+2 = O z k+j+3 . k+j+2 10 Sao cấp chính xác Khái niệm sao cấp chính xác được ra đưa khi nghiên cứu về tính ổn định của xấp xỉ Padé của hàm ez (Wanner, Hairer và Norsett 1978). Định nghĩa 1.6. Cho tập hợp A = {z ∈ C : |R (z)| > |ez |} = {z ∈ C : |q (z)| > 1} , ở đây q (z) = (1.17) R (z) . Tập A khi đó được gọi là sao cấp chính xác của R (z). ez Sao cấp chính xác không so sánh |R (z)| với 1 như khi xét tính ổn định, mà nó so sánh |R (z)| với nghiệm chính xác của phương trình thử |ez | = ex (z = x + iy) và hi vọng sẽ nhận được nhiều thông tin hơn. Chúng ta luôn giả sử rằng các hệ số của hàm R (z) là số thực và sao cấp chính xác đối xứng qua trục thực. Hơn nữa, với z = iy ta có |ez | = 1, khi đó tập A là phần bù của miền ổn định tuyệt đối S trên trục ảo. Hình 1.1 và Hình 1.2 thể hiện các tập sao cấp chính xác của xấp xỉ Padé. Bổ đề 1.1 (xem [3]). Nếu R (z) là một xấp xỉ cấp p của ez , tức là
z p+1 p+2 e − R (z) = Cz +O z với C 6= 0 thì khi z → 0, A có dáng điệu như một hình ngôi sao với p + 1 cánh sao có các góc quét bằng nhau và bằng 1.2 π . p+1 Xây dựng các phương pháp Runge-Kutta ẩn Phần này sẽ đề cập đến lớp các phương pháp Runge-Kutta ẩn sở hữu tính ổn định khá tốt. Việc xây dựng các phương pháp như vậy chủ yếu dựa vào các bộ điều kiện  s P 1   B (p) : bi cq−1 = q = 1, ..., p;  i  q  i=1   s P cqi q−1 i = 1, ..., s, q = 1, ..., η; C (η) : aij cj = (1.18) q  j=1   s  P bj  q
 bi cq−1 a = 1 − c j = 1, ..., s, q = 1, ..., ζ.  D (ζ) : ij i j i=1 q Điều kiện B (p) có nghĩa là công thức cầu phương (bi , ci ) có cấp p. Sự quan trọng của hai điều kiện còn lại được thể hiện ở Định lý 1.5 (xem [2] trang 208). 11 Hình 1.1: Tập sao cấp chính xác của xấp xỉ Padé với k = 1, j = 2. Hình 1.2: Tập sao cấp chính xác của xấp xỉ Padé với k = 2, j = 21. Định lý 1.5 (Butcher 1964). Nếu bộ hệ số bi , ci , aij của phương pháp RungeKutta thỏa mãn B (p) , C (η) , D (ζ) với p ≤ η + ζ + 1 và p ≤ 2η + 2 thì phương pháp có cấp p. 12 Các phương pháp Gauss Các phương pháp Gauss hay còn gọi "các phương pháp Kuntzmann-Butcher" là các phương pháp trùng khớp dựa trên cơ sở của công thức cầu phương Gauss với các hệ số c1 , ..., cs là nghiệm của đa thức trực giao Legendre bậc s ds s
s dxs x (x − 1) . Một số phương pháp Gauss được thể hiện trong Bảng 1.5. (a) p = 2 1 2 1 2 1 √ 1 3 − 2 √6 1 3 + 2 6 (b) p = 4 √ 1 3 − 4 6 1 4 1 2 1 4√ 1 3 + 4 6 1 2 Bảng 1.5: Các phương pháp Gauss cấp 2 và cấp 4. Định lý 1.6 (Butcher 1964, Ehle 1968). Phương pháp Gauss s nấc có cấp chính xác p = 2s. Hàm ổn định là xấp xỉ Padé (s, s) và phương pháp ổn định A. Chứng minh. Chi tiết xem [2] và [3]. Các phương pháp Radau IA và Radau IIA Butcher (1964) đã giới thiệu các phương pháp Runge-Kutta dựa trên cơ sở của các công thức cầu phương Radau và Lobatto. Ông ấy gọi chúng theo ba loại I, II hoặc III tùy thuộc vào c1 , c2 , ..., cs là nghiệm của I: II : III :
ds−1 xs (x − 1)s−1 , s−1 dx
ds−1 xs−1 (x − 1)s , s−1 dx
ds−2 xs−1 (x − 1)s−1 . s−2 dx Radau trái
Radau phải (Lobatto)
(1.19) (1.20) (1.21) Các trọng số bi được chọn để công thức cầu phương thỏa mãn B (s), trong đó B (2s − 1) là cơ sở Radau, B (2s − 2) là cơ sở Lobatto. Ehle (1969) đã theo đuổi ý tưởng của Butcher và xây dựng các phương pháp loại I, II và III với các tính chất ổn định vượt trội. Một cách độc lập, Axelsson (1969) đã tìm ra các phương pháp Radau IIA và chứng minh tính ổn định A của chúng. Phương pháp Radau IA s nấc thuộc loại I khi các hệ số aij (i, j = 1, ..., s) được định nghĩa bởi điều kiện D (s). Đây là cách duy nhất bởi vì các ci là phân biệt 13 (a) p = 1 0 (b) p = 3 1 1 0 2 3 1 4 1 4 1 4 (c) p = 5 1 − 4 5 12 3 4 √ −1 − 6 18 √ 88 + 7 6 360 √ 88 + 43 6 360 √ 16 + 6 36 1 9 1 9 1 9 1 9 0 √ 6− 6 10 √ 6+ 6 10 √ −1 + 6 18 √ 88 − 43 6 360 √ 88 − 7 6 360 √ 16 − 6 36 Bảng 1.6: Một số phương pháp Radau IA và bi 6= 0. Bảng 1.6 trình bày các phương pháp đầu tiên có đặc điểm này. Công thức loại II của Ehle thu được bằng cách áp dụng điều kiện C (s). Theo Định lý II.7.7 trong [2], các hệ số của phương pháp loại II được xây dựng dựa trên cơ sở là nghiệm của (1.20). Ta gọi chúng là phương pháp Radau IIA. Với s = 1 chúng ta thu được công thức Euler ẩn (xem Bảng 1.3). Các ví dụ về phương pháp Radau IIA được đưa ra trong Bảng 1.7. (a) p = 3 1 3 1 5 12 3 4 3 4 −1 12 1 4 1 4 √ 4− 6 10√ 4+ 6 10 1 (b) p = 5 √ 88 − 7 6 360 √ 269 + 169 6 1800 √ 16 − 6 36√ 16 − 6 36 √ 269 − 169 6 1800√ 88 + 7 6 360√ 16 + 6 36√ 16 + 6 36 √ −2 + 3 6 225 √ −2 − 3 6 225 1 9 1 9 Bảng 1.7: Các phương pháp Radau IIA cấp 3 và cấp 5. Định lý 1.7. Phương pháp Radau IA s nấc và phương pháp Radau IIA s nấc đều có cấp chính xác p = 2s − 1. Hàm xấp xỉ của chúng là xấp xỉ Padé (s − 1, s). Cả hai phương pháp đều ổn định A. Chứng minh. Chi tiết xem trong [3]. Các phương pháp Lobatto IIIA, IIIB và IIIC Với tất cả các công thức loại III thì ci là nghiệm của đa thức cho bởi (1.21) và các trọng số bi thỏa mãn B (2s − 2). Các hệ số aij được định nghĩa bởi C (s) cho ta các công thức Lobatto IIIA. Do đó nó là một phương pháp trùng khớp. Với các phương pháp Lobatto IIIB chúng ta áp dụng D (s). Cuối cùng, để có công thức Lobatto IIIC chúng ta đặt ai1 = b1 i = 1, 2, ..., s 14 (1.22) và xác định các aij còn lại bởi C (s − 1). Ehle (1969) đã giới thiệu hai lớp phương pháp đầu tiên, và trình bày các phương pháp IIIC với s ≤ 3. Định nghĩa chung của các phương pháp IIIC được đưa ra bởi Chipman (1971) và Axelsson (1972). Các ví dụ được đưa ra ở các Bảng 1.8 - 1.10. Định lý 1.8. Các phương pháp Lobatto IIIA, IIIB và IIIC có cấp chính xác p = 2s − 2. Hàm ổn định của các phương pháp Lobatto IIIA và Lobatto IIIB là xấp xỉ Padé (s − 1, s − 1). Với các phương pháp Lobatto IIIC, hàm ổn định là xấp xỉ Padé (s − 2, s). Và tất cả chúng đều ổn định A. Chứng minh. Chi tiết xem Phần IV.5 trong [3]. (a) p = 2 0 0 1 2 1 2 1 (b) p = 4 0 1 2 0 1 2 1 2 1 1 0 5 24 1 6 1 6 0 1 3 2 3 2 3 0 −1 24 1 6 1 6 Bảng 1.8: Các phương pháp Lobatto IIIA cấp 2 và cấp 4. (a) p = 2 0 1 1 2 1 2 1 2 (b) p = 4 0 0 0 1 2 1 2 1 1 6 1 6 1 6 1 6 −1 6 1 3 5 6 2 3 0 0 0 1 6 Bảng 1.9: Các phương pháp Lobatto IIIB cấp 2 và cấp 4. (a) p = 2 0 1 1 2 1 2 1 2 (b) p = 4 −1 2 1 2 1 2 0 1 2 1 1 6 1 6 1 6 1 6 −1 3 5 12 2 3 2 3 1 6 −1 12 1 6 1 6 Bảng 1.10: Các phương pháp Lobatto IIIC cấp 2 và cấp 4. Tóm tắt về việc xây dựng các phương pháp Runge-Kutta ẩn được thể hiện trong Bảng 1.11. 15 Phương pháp Bộ điều kiện đơn giản hóa Cấp chính xác Hàm ổn định Gauss B (2s) C (s) D (s) 2s (s, s) Padé Radau IA B (2s − 1) C (s − 1) D (s) 2s − 1 (s − 1, s) Padé Radau IIA B (2s − 1) C (s) D (s − 1) 2s − 1 (s − 1, s) Padé Lobatto IIIA B (2s − 2) C (s) D (s − 2) 2s − 2 (s − 1, s − 1) Padé Lobatto IIIB B (2s − 2) C (s − 2) D (s) 2s − 2 (s − 1, s − 1) Padé Lobatto IIIC B (2s − 2) C (s − 1) D (s − 1) 2s − 2 (s − 2, s) Padé Bảng 1.11: Thể hiện đầy đủ của các phương pháp Runge-Kutta ẩn. W -biến đổi Ta kí hiệu √ Pk (x) =  k k k 2k + 1 d x (x − 1)  dxk k! k X √ (−1)j+k = 2k + 1  k j  j+k j  xj j=0 (1.23) là đa thức Legendre đã được chuẩn hóa Z1 Pk2 (x) dx = 1. (1.24) 0 Các đa thức thỏa mãn Zx 1 P0 (t) dt = ξ1 P1 (x) + P0 (x) ; 2 (1.25a) 0 Zx Pk (t) dt = ξk+1 Pk+1 (x) − ξk Pk−1 (x) k = 1, 2... (1.25b) 0 với 1 ξk = √ 2 4k 2 − 1 (1.26) Định lý 1.9 (xem [3]). Cho W định nghĩa bởi wij = Pj−1 (ci ) i = 1, ..., s j = 1, ..., s, (1.27) và A là ma trận hệ số của phương pháp Gauss có cấp p = 2s thì   1/2 −ξ1 0 −ξ2  ξ1  W−1 AW =    ξ2 ... ... 16 ... 0 −ξs−1 ξs−1 0    := XG .   (1.28) Bổ đề 1.2 (xem [3]). Cho A là ma trận hệ số của một phương pháp Runge-Kutta ẩn và W là ma trận không suy biến với wij = Pj−1 (ci ) i = 1, ..., s j = 1, ..., η + 1. Khi đó điều kiện C (η) (với η ≤ s − 1) tương đương với điều kiện η cột đầu tiên của W−1 AW giống với η cột đầu tiên của XG trong (1.28). Bổ đề 1.3 (xem [3]). Cho W là ma trận không suy biến với wij = Pj−1 (ci ) i = 1, ..., s j = 1, ..., ζ + 1, và B = diag (b1 , ..., bs ) với bi 6= 0. Khi đó điều kiện D (ζ) (với ζ ≤ s − 1) tương

−1 đương với điều kiện ζ hàng đầu tiên của ma trận WT B A WT B giống với ζ hàng đầu tiên của XG trong (1.28). (Nếu B suy biến thì ta vẫn có (1.29) như bên dưới).   1/2 −ξ 1 ...   ξ1 0   ... ...  T  T −ξζ−1 (1.29) W BA =   W B.   ξ 0 −ξ ζ−1 ζ   ∗ ∗ ∗ ∗ ... ... ... ... ... ... ∗ ∗ Hai ma trận biến đổi ở trong các Bổ đề 1.2 và Bổ đề 1.3 ở trên có thể bằng nhau, tức là WT B = W−1 hoặc WT BW = I. (1.30) Bổ đề 1.4 (xem [3]). Với bất kì công thức cầu phương mà có cấp p ≥ 2s − 2 thì ma trận W = (Pj−1 (ci ))i,j=1,...,s (1.31) thỏa mãn (1.30). Định nghĩa 1.7. Với η , ζ là các số nguyên từ 0 đến s − 1. Ta nói rằng ma trận W cấp s × s thỏa mãn T (η, ζ) với công thức cầu phương có bộ hệ số (bi , ci )si=1 nếu   a) W là ma trận không suy biến    b) wij = Pj−1 (ci ) i = 1, ..., s, j = 1, ..., max (η, ζ) + 1 T (η, ζ)    I 0  T  c) W BW =  0 R 17 I là ma trận đơn vị cấp (ζ + 1) × (ζ + 1), R là ma trận tùy ý cấp (s − ζ − 1) × (s − ζ − 1) . Định lý 1.10 (xem [3]). Với W thỏa mãn T (η, ζ) cho công thức cầu phương (bi , ci )si=1 thì một phương pháp Runge-Kutta dựa trên (bi , ci )si=1 , với ma trận X = W−1 AW ta có a) η cột đầu tiên của X là của XG ⇔ C (η); b) ζ hàng đầu tiên của X là của XG ⇔ D (ζ). Bổ đề 1.5 (xem [3]). Nếu công thức cầu phương có các nút ci , tất cả các trọng số bi > 0 và nếu nó có cấp chính xác p thỏa mãn p ≥ 2η + 1, p ≥ 2ζ + 1 thì ma trận W = (pj−1 (ci ))i,j=1,...,s (1.32) có tính chất T (η, ζ) và thỏa mãn (1.30). Ở đây, pj (x) là đa thức trực giao bậc j với tích vô hướng s X hp, ri = bi p (ci ) r (ci ). i=1 Nhận xét 1.2 (xem [3]). Chúng ta có thể biểu diễn hàm ổn định của phương pháp Runge-Kutta ẩn dưới dạng ma trận Runge-Kutta chuyển vị X = W−1 AW. Từ (b) và (c) của tính chất T (η, ζ), ta suy ra We1 = 1 , WT B11 = e1 , e1 = (1, 0, ..., 0)T . (1.33) Do đó các công thức (1.10) và (1.11)trở thành R (z) = 1 + zeT1 (I − zX)−1 e1 , (1.34)
T R (z) = 1.3 det I − zX + ze1 e1 . det (I − zX) (1.35) Áp dụng các phương pháp Runge-Kutta giải bài toán cương Để minh họa cho các phương pháp Runge-Kutta được xây dựng trong phần trên, ta áp dụng các phương pháp Runge-Kutta giải một số bài toán cương. Trước hết ta tìm hiểu xem bài toán cương là gì? Có nhiều cách định nghĩa khác nhau về bài toán cương. Tuy nhiên, quan điểm có tính thực tiễn nhất 18 cũng là quan điểm đầu tiên trong lịch sử về bài toán cương là do Curtiss và Hirschfelder đưa ra. Theo đó, bài toán cương là những bài toán mà các phương pháp Runge-Kutta hiển không giải quyết được. Đối với các bài toán cương, các phương pháp Runge-Kutta ẩn giải quyết tốt hơn rất nhiều so với các phương pháp Runge-Kutta hiển. Định nghĩa 1.8. Bài toán giá trị ban đầu được gọi là bài toán cương nếu điều kiện ổn định tuyệt đối xác định bước đi h nhỏ hơn nhiều so với yêu cầu về độ chính xác. Ngoài ra, ta có thể nhận dạng bài toán cương trong một vài trường hợp đặc biệt. Với bài toán (1.6) ta có định nghĩa sau Định nghĩa 1.9. Đối với bài toán y 0 = Ay, A ∈ Cm×m nếu Re (λi ) ≤ 0 (với λi là max |Reλi | các giá trị riêng của A) với i = 1, ..., m và i min |Reλi |  1 thì bài toán gọi là bài i toán cương. Bài toán 1.1. Bài toán Van der Pol  0  y1 = y2

y20 = 1 − y12 y2 − y1 /ε  0 ≤ x ≤ 3. (1.36) y1 (0) = 2, y2 (0) = 0 Ta xét bài toán với ε = 10−3 . Ta lập trình chương trình thử nghiệm số để giải bài toán Van der Pol bằng các phương pháp Gauss cấp 2 và cấp 4 (xem Bảng 1.5) trong môi trường Matlab. Sau đó, ta so sánh kết quả lời giải số của chương trình trên với lời giải số của chương trình ode23s. Chương trình ode23s được xây dựng dựa trên cặp công thức Rosenbrock cấp 2 và cấp 3 (xem [6]). Kết quả các lời giải số của y1 thể hiện trong Hình 1.3. Sau 2000 bước với bước đi đều h = 0.0015, phương pháp Gauss cấp 4 cho ta lời giải khá chính xác, còn lời giải của phương pháp Gauss cấp 2 có nhiều sai số. Trong khi đó, nếu ta sử dụng chương trình ode23s thì chỉ cần 450 bước với hmin ≈ 1.4629 × 10−6 , hmax ≈ 0.0431 ta đã có lời giải khá chính xác. Bài toán 1.2. Bài toán cương Robertson  0 y1 = −0.04y1 + 104 y2 y3    y 0 = 0.04y − 104 y y − 3.107 y 2 2  y0 =   3 1 2 3 2 3.107 y22 y1 (0) = 1; y2 (0) = 0; y3 (0) = 0, 0 ≤ x ≤ 40. 19 (1.37)