Tài liệu Phương pháp số giải phương trình vi phân có chậm

1 ĐẠI HỌC QUỐC GIA HÀ NỘI TRƯỜNG ĐẠI HỌC KHOA HỌC TỰ NHIÊN ĐỖ THỊ THÚY NGỌC PHƯƠNG PHÁP SỐ GIẢI PHƯƠNG TRÌNH VI PHÂN CÓ CHẬM Chuyên ngành : TOÁN HỌC TÍNH TOÁN Mã số : 60 46 30 Người hướng dẫn : PGS. TS. VŨ HOÀNG LINH HÀ NỘI - 2012 2 Mục lục Lời nói đầu 4 1 Giới thiệu 1.1 Một vài ví dụ so sánh phương trình vi phân có chậm và phương trình vi phân thường . . . . . . . . . . . . . . . . . . . . . . . . . . 1.2 Phương pháp số giải phương trình vi phân thường . . . . . . . . 1.2.1 Các khái niệm cơ bản . . . . . . . . . . . . . . . . . . . . . 1.2.2 Một số phương pháp số tiêu biểu giải phương trình vi phân thường . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.3 Nghiệm số của phương trình vi phân có chậm: Phương pháp cho phương trình vi phân thường liệu có đủ hay không? . . . . . . . . 1.3.1 Sự thất bại về cấp chính xác của phương pháp . . . . . . . 1.3.2 Sự thất bại về tính ổn định . . . . . . . . . . . . . . . . . 1.3.3 Một phương pháp tốt cho các PTVPCC . . . . . . . . . . 6 2 Sự tồn tại và tính chính qui của nghiệm của phương trình vi phân có chậm 2.1 Vị trí của các điểm gián đoạn và sự trơn dần của nghiệm . . . . 2.1.1 Các điểm gián đoạn gốc và thứ cấp . . . . . . . . . . . . . 2.1.2 Chậm triệt tiêu và không triệt tiêu . . . . . . . . . . . . . 2.1.3 Chậm bị chặn và không bị chặn . . . . . . . . . . . . . . . 2.2 Sự tồn tại và tính duy nhất nghiệm . . . . . . . . . . . . . . . . . . . . 7 . . . 12 . . . 12 . . . 14 . . . . . . . . . . . . . . . . . . . . . . 17 18 20 25 . . . . . 27 27 28 31 32 35 3 Các phương pháp cho phương trình vi phân có chậm 37 3.1 Hướng tiếp cận đầu tiên . . . . . . . . . . . . . . . . . . . . . . . . . . . 38 3.2 Các kết quả sơ bộ về phương pháp số giải PTVP với đầu ra liên tục . 40 3.3 Hướng tiếp cận thông thường thông qua phương pháp số giải PTVP với đầu ra liên tục . . . . . . . . . . . . . . . . . . . . . . . . . . 47 4 Sự hội tụ của phương pháp số giải PTVP với đầu ra liên tục 50 3 4.1 PTVPCC với chậm hằng số hoặc chậm phụ thuộc thời gian không triệt tiêu . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51 4.2 PTVPCC với chậm phụ thuộc thời gian tuỳ ý . . . . . . . . . . . . . . 59 Kết luận 67 Tài liệu tham khảo 69 Phụ lục 70 4 LỜI NÓI ĐẦU Ta đã biết phương trình vi phân thường (PTVPT) là phương trình có dạng: y 0 (t) = f (t, y(t)), t ≥ t0 . Tuy nhiên, trong nhiều vấn đề, biến y 0 còn phụ thuộc vào các giá trị trong quá khứ của biến y . Khi đó PTVPT trên biến đổi thành phương trình sau, gọi là phương trình vi phân có chậm (PTVPCC): y 0 (t) = f (t, y(t − τ1 ), ..., y(t − τn )), t ≥ t0 , với τi = τi (t, y(t)) ≥ 0 ∀t ≥ t0 , i = 1, ..., n, được gọi là các chậm. Việc nghiên cứu về mặt lí thuyết cũng như phương pháp số giải PTVPT đã thu hút được sự quan tâm của các nhà toán học trong một thời gian dài với rất nhiều những kết quả quan trọng. Ngược lại, việc nghiên cứu đối với PTVPCC mới được quan tâm nhiều trong thời gian từ những thập niên cuối thế kỉ 20 trở lại đây. Sự quan tâm của các nhà toán học ứng dụng dành cho các phương pháp số giải PTVPCC ngày càng gia tăng, thể hiện qua số lượng ngày càng nhiều các sách chuyên khảo, bài báo và công trình nghiên cứu được công bố và đăng tải trên các tạp chí toán học uy tín. Mục tiêu của luận văn này là giới thiệu các phương pháp số giải PTVPCC, đặc biệt là các phương pháp số giải PTVP với đầu ra liên tục áp dụng giải bài toán giá trị ban đầu cho PTVP có chậm hằng số hoặc chậm chỉ phụ thuộc thời gian. Luận văn gồm 4 chương: Chương 1 sẽ giới thiệu PTVPCC thông qua việc so sánh với PTVPT. Những khác biệt quan trọng nhất về tính chất định tính cũng như về khía cạnh giải số sẽ được đề cập đến trong phần đầu của chương. Phần thứ hai của chương giới thiệu một cách vắn tắt các khái niệm cơ bản về phương pháp số giải PTVPT, từ đó phần thứ ba sẽ chỉ ra lí thuyết cho PTVPT là không đủ khi áp dụng cho PTVPCC. 5 Chương 2 trình bày sự tồn tại và tính chính qui nghiệm của PTVPCC. Đặc biệt, ta sẽ tập trung sự quan tâm vào tính chất và vị trí của các điểm gián đoạn của đạo hàm, nếu có, và sự lan truyền của chúng dọc theo khoảng tích phân dưới các giả thiết khác nhau trên các chậm. Chương 3 phác hoạ một cách vắn tắt một vài hướng tiếp cận khác nhau đã được sử dụng để giải số các PTVPCC, tập trung vào hướng tiếp cận thông thường thông qua các phương pháp số giải PTVP với đầu ra liên tục. Chương 4 trình bày sự hội tụ của các phương pháp số giải PTVP với đầu ra liên tục đã được giới thiệu trong chương 3. Đặc biệt, ta chứng minh tính đặt chỉnh (well-posedness) của phương pháp và phân tích cấp hội tụ của nó cho PTVPCC và PTVPCC trung tính có chậm hằng số hoặc chậm phụ thuộc thời gian. Luận văn được hoàn thành trên cơ sở tham khảo hai tài liệu chính: 1. Uri M.Ascher, Linda R.Petzold, Computer Methods for Ordinary Equations and Differential-Algebraic Equations, SIAM, 1998. 2. Alfredo Bellen, Marino Zennaro, Numerical Methods for Delay Differential Equations, Oxford University Press, 2003. Qua đây, tác giả xin chân thành cảm ơn PGS.TS Vũ Hoàng Linh, người thầy đã tận tâm giảng dạy và hướng dẫn tác giả hoàn thành bản luận văn này. Xin được gửi lời cảm ơn chân thành đến trường Đại học Khoa học Tự nhiên - Đại học Quốc gia Hà Nội nơi tác giả đã hoàn thành chương trình cao học dưới sự giảng dạy và hướng dẫn nhiệt tình của các thầy các cô. Xin được cảm ơn tập thể lãnh đạo, chuyên viên phòng Giáo dục Trung học, Sở Giáo dục và Đào tạo Ninh Bình đã thông cảm, động viên và tạo mọi điều kiện thuận lợi về thời gian, công việc để tác giả hoàn thành khoá học và luận văn. Cuối cùng, xin được cảm ơn bè bạn và gia đình đã hỗ trợ, động viên và chia sẻ những khó khăn với tác giả trong suốt thời gian học tập vừa qua. Do thời gian và trình độ còn hạn chế, chắc chắn bản luận văn không thể tránh khỏi những thiếu sót, tác giả rất mong nhận được sự chỉ bảo tận tình của các thầy cô và bạn bè đồng nghiệp, tác giả xin chân thành cảm ơn! Hà Nội, 15.11.2012 Học viên Đỗ Thị Thuý Ngọc 6 CHƯƠNG 1 GIỚI THIỆU Nhiều vấn đề thực tế trong vật lí, kĩ thuật, sinh học, y học, kinh tế . . . có thể được mô hình hóa bằng một bài toán giá trị ban đầu, hoặc còn được gọi là bài toán Cauchy, cho các PTVPT có dạng  y 0 (t) = g(t, y(t)), t ≥ t0 , (1.1) y(t ) = y , 0 0 trong đó hàm y(t), được gọi là biến trạng thái, biểu diễn một đại lượng nào đó tham gia vào quá trình. Tuy nhiên, để làm cho mô hình phù hợp hơn với các hiện tượng thực tế, đôi khi ta cần biến đổi vế phải của (1.1) để thể hiện sự phụ thuộc của biến y 0 vào các giá trị trong quá khứ của biến trạng thái y . Dạng tổng quát nhất của các mô hình như thế được cho bởi PTVPCC y 0 (t) = f (t, yt ), t ≥ t0 , trong đó yt = y(t + θ), θ ∈ [−r, 0], là một hàm thuộc vào không gian Banach C = C 0 ([−r, 0], Rd ) các hàm số liên tục ánh xạ khoảng [−r, 0] vào Rd , và f : Ω → Rd là một hàm đã cho ánh xạ tập hợp Ω ⊂ R × C vào Rd . Bài toán giá trị ban đầu bây giờ là  y 0 (t) = f (t, yt ), t ≥ t0 , (1.2) y = y(t + θ) = Φ(θ), t0 0 trong đó Φ(θ) ∈ C biểu diễn điểm khởi tạo hoặc dữ liệu khởi tạo. 7 1.1. Một vài ví dụ so sánh phương trình vi phân có chậm và phương trình vi phân thường Trong tài liệu này, bài toán giá trị ban đầu (1.2) sẽ được mô tả theo một cách thức thân thiện hơn như sau  y 0 (t) = f (t, y(t − τ1 ), ...y(t − τn )), t ≥ t0 , (1.3) y(t) = φ(t), t ≤ t . 0 Tùy theo độ phức tạp của hiện tượng, các chậm τi luôn luôn là không âm, có thể chỉ là các hằng số (trường hợp chậm hằng số), hoặc các hàm số của t, τi = τi (t) (trường hợp chậm biến thiên hoặc phụ thuộc thời gian) hoặc thậm chí là các hàm số của t và chính y , τi = τi (t, y(t)) (trường hợp chậm phụ thuộc trạng thái). Luận văn này sẽ tập trung sự nghiên cứu vào hai trường hợp: chậm phụ thuộc hằng số và chậm phụ thuộc thời gian. Để đơn giản hóa về mặt kí hiệu, hàm φ(t) được hiểu là được định nghĩa trong [ρ, t0 ], trong đó   ρ = min 1≤i≤n min(t − τi ) . t≥t0 Một trường hợp khá phổ biến và thú vị là khi n = 2 và τ1 ≡ 0, khi đó (1.3) có dạng như sau  y 0 (t) = f (t, y(t), y(t − τ )), t ≥ t0 , (1.4) y(t) = φ(t), t ≤ t . 0 Vì với t ≥ t0 nào đó có thể xảy ra t − τ < t0 , sự khác nhau đầu tiên giữa các phương trình dạng (1.1) và (1.4) đó là nghiệm của (1.4) thường phụ thuộc vào hàm khởi tạo φ(t) hơn là phụ thuộc vào giá trị khởi tạo y0 như đối với (1.1). Nói chung, đạo hàm bên phải y 0 (t+ 0 ), đó là f (t0 , φ(t0 ), φ(t0 − τ )), không bằng đạo hàm − trái y 0 (t0 ) và do đó nghiệm y không được liên kết trơn với hàm khởi tạo φ(t) tại điểm t0 , ở đó chỉ có tính chất C 0 -liên tục là có thể được đảm bảo. Hơn nữa, tính không liên tục của đạo hàm sẽ lan truyền từ điểm khởi tạo t0 theo khoảng tích phân và tạo ra các điểm gián đoạn tiếp theo mà tại đó, nghiệm càng ngày càng trơn hơn. Như một hệ quả, thậm chí nếu các hàm f (t, y, x), τ (t, y) và φ(t) trong (1.4) là C ∞ -liên tục thì nói chung y(t) đơn giản chỉ là C 1 -liên tục trong [t0 , tf ]. Ví dụ 1.1.1. Xét phương trình  y 0 (t) = −y(t − 1), t ≥ 0, y(t) = 1, t ≤ 0. (1.5) 8 Nghiệm của phương trình được miêu tả trong Hình 1.1. Vì y 0 (0− ) = 0 và y 0 (0+ ) = −y(−1) = −1, hàm đạo hàm y 0 (t) có một điểm gián đoạn tại t = 0. Đạo hàm cấp hai y 00 (t) được cho bởi y 00 (t) = −y 0 (t − 1), và do đó nó có một điểm gián đoạn tại t = 1. Đạo hàm cấp ba được cho bởi y 000 (t) = −y 00 (t − 1) = y 0 (t − 2) và do đó nó có một điểm gián đoạn tại t = 2, và tương tự như thế tại các điểm là các bội số của chậm t = 3, 4, . . . Hình 1.1: Nghiệm của (1.5) Trong các mô hình tổng quát hơn, đạo hàm y 0 (t) có thể phụ thuộc vào y và chính y 0 tại một giá trị quá khứ t − τ nào đó. Trong trường hợp này, (1.4) thay đổi thành dạng  y 0 (t) = f (t, y(t), y(t − τ ), y 0 (t − τ )), t ≥ t0 , (1.6) y(t) = φ(t), t ≤ t , 0 trong đó hàm φ(t) được giả thiết ít nhất là C 1 -liên tục. Phương trình (1.6) được gọi là một PTVPCC trung tính (delay differential equation of neutral type). Như trước đã nói, hàm khởi tạo φ(t) không liên kết trơn với nghiệm y(t) tại t0 , nơi chỉ có duy nhất tính liên tục được đảm bảo. Điểm gián đoạn này sẽ lan truyền thành một tập các điểm gián đoạn mà ở đó nghiệm, không giống trường hợp không trung tính, chỉ thuộc duy nhất lớp C 0 . Do đó, trừ phi điều kiện nối 0 φ0 (t− 0 ) = f (t0 , φ(t0 ), φ(t0 − τ ), φ (t0 − τ )) được thỏa mãn, còn không thì nghiệm của (1.6) phải được hiểu theo nghĩa tổng quát “hầu khắp nơi”. 9 Ví dụ 1.1.2. Xét phương trình  y 0 (t) = −y 0 (t − 1), t ≥ 0, y(t) = t, t ≤ 0, (1.7) có nghiệm được vẽ trên Hình 1.2. Vì y 0 (0− ) = 1 và y 0 (0+ ) = −y 0 (−1) = −1, đạo hàm y 0 (t) có một điểm gián đoạn tại t = 0. Hơn nữa, vì y 0 (t) = −y 0 (t − 1) với mọi t ≥ 0, đạo hàm y 0 (t) không liên tục tại t = 1 cũng như tại t = 2, 3, . . . là các bội số của t = 1. Hình 1.2: Nghiệm của (1.7). Ví dụ sau chỉ ra rằng, trong khi các nghiệm bị chặn của các PTVPT có thể dao động chỉ khi hệ thống có ít nhất hai thành phần và có thể dao động hỗn loạn chỉ khi hệ thống có ít nhất ba thành phần (Định lí Poincaré-Bendixon), thì các nghiệm của PTVPCC có thể có tính chất dao động và thậm chí dao động hỗn loạn trong trường hợp vô hướng. Ví dụ 1.1.3. Xét phương trình logistic có chậm sau y 0 (t) = ay(t)(1 − y(t − 1)), (1.8) phương trình này mô hình hóa sự thay đổi của dân số, là cải tiến của mô hình Verhulst-Pearl y 0 (t) = ay(t)(1 − y(t)). Trong khi các nghiệm của phương trình Verhulst-Pearl là đơn điệu, các nghiệm dương của (1.8) là đơn điệu với a ∈ (0, 1/e), dao động với a ∈ [1/e, π/2) và xấp xỉ với các quỹ đạo tuần hoàn với a > π/2 (Hình 1.3 và 1.4). (Xem [6]) Cuối cùng, ta thấy rằng sự có mặt của thành phần chậm có thể thay đổi mạnh mẽ tính chất định tính của nghiệm bằng cách tác động đến sự ổn định của mô hình. 10 Hình 1.3: Nghiệm của (1.8) với y(t) = 0.1 khi t ≤ 0, a = 1.4 và 0.3. Hình 1.4: Nghiệm của (1.8) với a = 1.7 trong mặt phẳng pha. Ví dụ 1.1.4. Xét phương trình vô hướng tuyến tính  y 0 (t) = λy(t) + µy(t − 1), t ≥ 0, y(t) = −t + 1, t ≤ 0, (1.9) với các hệ số λ, µ là các hằng số thực. Ta biết rằng, với µ = 0, phương trình (1.9) có dạng  y 0 (t) = λy(t), t ≥ 0, (1.10) y(0) = 1, Nghiệm của phương trình này tiệm cận tới 0 với mọi λ âm và bùng nổ với λ dương bất kì. Hơn nữa trong trường hợp λ âm, nghiệm còn bị chặn bởi giá 11 trị khởi tạo 1. Mặt khác, với µ 6= 0, thành phần chậm µy(t − 1) trong (1.9) tác động như một thành phần cưỡng bức và các tính chất đã nêu của nghiệm có thể không còn được thoả mãn. Đặc biệt, với mọi µ > 0, tồn tại λ < 0 sao cho nghiệm không tiệm cận tới 0 và tồn tại giá trị λ < 0 khác sao cho nghiệm tiệm cận tới 0 nhưng không bị chặn bởi giá trị khởi tạo y(0) = 1. Các tình huống này được minh họa trong Hình 1.5 với µ = 4, λ = −3.5 và λ = −5. Cũng vậy, với λ = 0.5 và µ = −1, thành phần chậm −y(t − 1) tác động như là một thành phần ổn định hoá (stabilizer) của mô hình có nghiệm ổn định bất chấp tính dương của λ (xem Hình 1.6). Hình 1.5: Nghiệm ổn định và không ổn định của (1.9) với λ < 0. Hình 1.6: Nghiệm ổn định của (1.9) với λ > 0. Ví dụ 1.1.5. PTVPCC trung tính sau là một ví dụ cho thấy chậm với giá trị nhỏ có thể tạo ra một ảnh hưởng rộng lớn y 0 (t) = −1.5y 0 (t − τ ) + λy(t), λ < 0, (1.11) 12 Nghiệm của phương trình này ổn định tiệm cận với τ = 0 và không ổn định với mọi τ > 0. Trong trường hợp này thành phần chậm tác động như là một thành phần làm mất ổn định (destabilizer). 1.2. Phương pháp số giải phương trình vi phân thường Để có cơ sở theo dõi được điểm khác biệt giữa các phương pháp số giải PTVPT và PTVPCC, trong phần này ta sẽ nhắc lại các khái niệm cơ bản liên quan đến PTVPT, một số phương pháp số tiêu biểu giải PTVPT cũng như tính chất của các phương pháp đó. (Xem [3]) 1.2.1. Các khái niệm cơ bản Xét bài toán giá trị đầu cho PTVPT tổng quát  y 0 = f (t, y), 0 ≤ t ≤ b, y(0) = c, (1.12) trong đó b > 0, c ∈ Rn cho trước, f ∈ C([0, b] × Rn → Rn ). Định nghĩa 1.2.1. Nghiệm y(t) được gọi là ổn định nếu với mọi ε > 0 cho trước, tồn tại một δ > 0 sao cho với mọi nghiệm ŷ(t) khác y(t) và thoả mãn |y(0) − ŷ(0)| ≤ δ ta cũng có |y(t) − ŷ(t)| ≤ ε, ∀t ≥ 0. Định nghĩa 1.2.2. Nếu có thêm điều kiện |y(t) − ŷ(t)| → 0 khi t → ∞ thì y(t) được gọi là ổn định tiệm cận. Chọn lưới điểm ∆ = {t0 = 0, t1 , . . . , tN = b} và đặt hn = tn − tn−1 với n = 1, . . . , N . Ta xét một phương pháp số đơn giản giải (1.12), đó là phương pháp Euler hiển, có dạng yn = yn−1 + hn f (tn−1 , yn−1 ). Ta viết lại công thức trên thành yn − yn−1 − f (tn−1 , yn−1 ) = 0. hn 13 Cho u là một hàm bất kì định nghĩa trên lưới điểm ∆. Xét toán tử sai phân Nh u(tn ) ≡ u(tn ) − u(tn−1 ) − f (tn−1 , u(tn−1 )) hn với n = 1, . . . , N . Toán tử sai phân này thay đổi tùy theo phương pháp. Xét yh là một hàm lưới nhận giá trị yn tại mỗi điểm tn , n = 0, 1, . . . , N . Khi đó phương pháp số cho bởi phương trình toán tử Nh yh (tn ) = 0, n = 0, ..., N (với y0 = c). Định nghĩa 1.2.3. Gọi y là nghiệm chính xác của phương trình. Sai số chặt cụt địa phương là đại lượng được định nghĩa bởi dn = Nh y(tn ). Định nghĩa 1.2.4. Phương pháp số được gọi là chính xác cấp p nếu Nh y(tn ) = O (hpn ) với p là một số nguyên dương. Đặt h = max hn . 1≤n≤N Định nghĩa 1.2.5. Phương pháp số được gọi là hội tụ cấp p nếu sai số toàn cục en , trong đó en = yn − y(tn ), e0 = 0 thoả mãn en = O(hp ) với n = 1, 2, . . . , N . Định nghĩa 1.2.6. Phương pháp số được gọi là ổn định – 0 nếu có các hằng số h0 và K sao cho với các hàm lưới xh và zh bất kì với h ≤ h0 ta luôn có   |xn − zn | ≤ K |x0 − z0 | + max |Nh xh (tj ) − Nh zh (tj )| , 1 ≤ n ≤ N. 1≤j≤N Định lý 1.2.1. Nếu phương pháp là chính xác cấp p và ổn định – 0 thì nó hội tụ cấp p |en | ≤ K max |dj | = O (hp ) . j Xét phương trình thử y 0 = λy với λ ∈ C, Reλ < 0 và y0 , y1 , ..., yN là một lời giải số. 14 Định nghĩa 1.2.7. Điều kiện sau được gọi là điều kiện ổn định tuyệt đối |yn | ≤ |yn−1 | , n = 1, 2, ..., N. Định nghĩa 1.2.8. Miền ổn định tuyệt đối của một phương pháp số là miền trong mặt phẳng phức z sao cho khi áp dụng phương pháp cho phương trình thử, với z = λh nằm trong miền này, ta được một nghiệm xấp xỉ thoả mãn điều kiện ổn định tuyệt đối. Định nghĩa 1.2.9. Một phương pháp số được gọi là ổn định – A nều miền ổn định tuyệt đối của nó chứa toàn bộ nửa bên trái mặt phẳng phức z = λh. 1.2.2. Một số phương pháp số tiêu biểu giải phương trình vi phân thường Trong phần này ta xét hai loại phương pháp tiêu biểu để giải (1.12), đó là các phương pháp một bước và các phương pháp đa bước. Một lớp phương pháp một bước quan trọng để giải (1.12) là các phương pháp Runge – Kutta (R - K). Sau đây là một số phương pháp R - K đơn giản nhất. (a) Phương pháp Euler hiển yn = yn−1 + hn f (tn−1 , yn−1 ). Người ta chứng minh được rằng phương pháp Euler hiển chính xác cấp 1, ổn định – 0 và do đó hội tụ cấp 1. Tuy nhiên phương pháp Euler hiển không ổn định – A. (b) Phương pháp Euler ẩn yn = yn−1 + hn f (tn , yn ). Cũng như phương pháp Euler hiển, người ta chứng minh được rằng phương pháp Euler ẩn chính xác cấp 1, ổn định – 0 và do đó hội tụ cấp 1. Tuy nhiên, không giống phương pháp Euler hiển, phương pháp Euler ẩn là ổn định – A. (c) Phương pháp trung điểm yn = yn−1 + hn f (tn−1/2 , yn +y2 n−1 ). (d) Phương pháp hình thang yn = yn−1 + h2n [f (tn−1 , yn−1 ) + f (tn , yn )] . Người ta chứng minh được rằng phương pháp trung điểm và phương pháp hình thang đều chính xác cấp 2, ổn định – 0 và do đó hội tụ cấp 2. Ngoài ra, cả hai phương pháp đều ổn định – A. Nói chung, một phương pháp R – K s - nấc giải PTVPT y 0 = f (t, y) có thể được viết dưới dạng Yi = yn−1 + h yn = yn−1 + h s P j=1 s P i=1 aij f (tn−1 + cj h, Yj ), i = 1, 2, ..., s, bi f (tn−1 + ci h, Yi ), 15 trong đó  s P   aij , i = 1, ..., s,  ci = j=1 s P   bi = 1.  i=1 Nếu aij = 0 ∀j > i thì ta có một công thức R – K hiển. Ngược lại ta có một công thức R – K ẩn. Ta thể hiện phương pháp một cách thuận tiện bằng cách sử dụng bảng Butcher như sau c1 c2 a11 a12 . . . a1s a21 a22 . . . a2s cs as1 as2 . . . ass b1 b2 . . . b3 .. . .. . .. . ... .. . Chẳng hạn, ta có bảng Butcher của (a) Phương pháp Euler hiển 0 0 1 1 1 1 (b) Phương pháp Euler ẩn (c) Phương pháp trung điểm 1/2 1/2 1 (d) Phương pháp hình thang 0 1 0 0 1/2 1/2 1/2 1/2 Ngoài các phương pháp một bước, người ta còn xây dựng các phương pháp đa bước để giải bài toán giá trị ban đầu cho các PTVPT. Một phương pháp k bước tuyến tính tổng quát có dạng k X j=0 αj yn−j = h k X βj fn−j , j=0 trong đó các hằng số αj , βj , j = 0, . . . , k là các hệ số của phương pháp, α0 6= 0; |αk | + |βk | = 6 0 và tn−j = tn − jh, yn−j ≈ y(tn−j ), fn−j = f (tn−j , yn−j ) ≈ f (tn−j , y(tn−j )), j = 0, ..., k. 16 Nếu β0 = 0, ta có một công thức đa bước hiển. Ngược lại, ta có một công thức đa bước ẩn. Hai họ phương pháp đa bước phổ biến nhất là các phương pháp Adams và các phương pháp BDF. Các phương pháp Adams hiển, còn được gọi là các phương pháp AdamsBashforth, là phổ biến nhất trong các phương pháp đa bước hiển và có dạng như sau k X βj fn−j yn = yn−1 + h j=1 trong đó βj = (−1) j−1  k−1 P i=j−1 1 iR γi = (−1) 0  −s i i j−1  γi ,  ds. Các phương pháp Adams ẩn, còn được gọi là các phương pháp AdamsMoulton, có dạng k X yn = yn−1 + h βj fn−j . j=0 Cấp của một phương pháp Adams-Moulton k -bước là p = k + 1. Các phương pháp Adams-Moulton có các hằng số sai số nhỏ hơn, sử dụng ít hơn một bước và có miền ổn định rộng hơn nhiều so với các phương pháp Adams-Bashforth có cùng cấp. Một phương pháp k -bước BDF với cấp p = k có dạng sau k X 1 i=1 i ∇i yn =hf (tn , yn ) , trong đó ∇0 f (tn , yn ) = f (tn , yn ) , ∇i f (tn , yn ) = ∇i−1 f (tn , yn ) − ∇i−1 f (tn−1 , yn−1 ) . Phương pháp BDF có tính stiff decay (giảm nhanh), thích hợp giải bài toán cương, xem [3]. 17 1.3. Nghiệm số của phương trình vi phân có chậm: Phương pháp cho phương trình vi phân thường liệu có đủ hay không? Để minh hoạ một vài đặc trưng cơ bản của các phương pháp số giải PTVPCC và sự khác nhau của chúng so với các phương pháp số giải PTVPT, ta xem xét PTVPCC hằng số sau  y 0 (t) = f (t, y(t), y(t − 1)), t ≥ 0, (1.13) y(t) = φ(t), t ≤ 0. Hướng tiếp cận việc giải số (1.13) một cách tự nhiên nhất, nhưng không phải là duy nhất, là sử dụng các bước tích phân nhỏ hơn hoặc bằng chậm τ = 1 và tích phân từng bước các PTVPT đạt được từ (1.13) bằng cách thay thành phần chậm y(t − 1) bằng một hàm η(t − 1) cho trước, tuỳ theo giá trị của t mà bằng hàm khởi tạo φ(t−1) hoặc bằng một mở rộng liên tục của nghiệm tìm được trước đó bởi chính phương pháp đang sử dụng. Do đó, tại bước thứ (n + 1), phương trình được giải sẽ là  0 ωn+1 (t) = f (t, ωn+1 (t), x(t − 1)), tn ≤ t ≤ tn+1 , (1.14) ωn+1 (tn ) = yn , trong đó  x(s) = φ(s), s ≤ 0, η(s), 0 ≤ s ≤ tn . Công thức tích phân cho ta giá trị yn+1 và nghiệm xấp xỉ η của (1.13) liên tục trong [tn , tn+1 ] sao cho η(tn+1 ) = yn+1 . Một tính chất riêng của hướng tiếp cận này là trong khi phương pháp số giải PTVPT chỉ cung cấp các giá trị xấp xỉ của nghiệm tại các điểm nút thì việc cài đặt phương pháp số giải (1.14) có thể yêu cầu các hiểu biết về nghiệm xấp xỉ η(t) tại một vài điểm t − 1 có thể khác các điểm nút. Do đó, nói chung, các phương pháp số giải PTVPCC sẽ dựa trên cơ sở các mở rộng liên tục của các phương pháp số giải PTVPT. Điều này có thể thực hiện bằng một phép nội suy hậu nghiệm (a posteriori interpolation) của các giá trị yn được cung cấp bởi phương pháp PTVPT rời rạc cơ sở hoặc, tốt hơn, bằng phương pháp số giải PTVP với đầu ra liên tục, đó là các phương pháp cung cấp từng bước một xấp xỉ liên tục của nghiệm (xem Hình 1.7). Như chúng ta sẽ thấy, sự thành công của phương pháp số giải PTVP thường với đầu ra liên tục đạt được về mặt sự chính xác và 18 Hình 1.7: Nghiệm xấp xỉ của (1.13) đạt được bằng phương pháp số giải PTVP với đầu ra liên tục. tính ổn định phụ thuộc vào sự lựa chọn đặc biệt phương pháp rời rạc cũng như là mở rộng liên tục. Ta cũng đã chỉ ra rằng sự xuất hiện của thành phần có chậm có thể thay đổi một cách mạnh mẽ một vài tính chất về tính bị chặn và ổn định, và nói chung, về sự biến đổi của mô hình PTVPT ban đầu. Bây giờ, ta muốn minh hoạ, qua các ví dụ, rằng trong khi cài đặt, một vài tính chất về sự chính xác và tính ổn định được mong đợi của phương pháp PTVPT cơ sở có thể bị phá huỷ khi phương pháp được áp dụng cho một PTVPCC. 1.3.1. Sự thất bại về cấp chính xác của phương pháp Để minh hoạ cho khả năng có thể mất tính chính xác, xét lớp các phương trình tuyến tính hệ số hằng số sau  0 y (t) = ay(t) − π2 ea y(t − 1), t ≥ 0, (1.15) π at y(t) = φ(t) = e sin( 2 t), t ≤ 0, Nghiệm của phương trình, y(t) = eat sin( π2 t), thuộc lớp C ∞ trong [−1, +∞). Theo (1.14), với n = 0, 1, . . . , ta sẽ giải PTVPT sau  0 wn+1 (t) = awn+1 (t) − π2 ea x(t − 1), tn ≤ t ≤ tn+1 , (1.16) wn+1 (tn ) = yn , trong đó  x(s) = φ(s) = eas sin π2 s , s ≤ 0, η(s), 0 ≤ s ≤ tn
Một lớp phương pháp tốt để tích phân (1.16) được cho bởi phương pháp trùng khớp tại ν điểm Gauss, có thể được xem là các phương pháp Runge-Kutta 19 ν -nấc cấp 2ν . Vì chúng là các phương pháp trùng khớp dựa trên sự xấp xỉ đa thức từng khúc bậc ν , nên chúng cũng cung cấp một mở rộng liên tục η(t) có cấp chính xác đều ν + 1. Do đó, phương pháp trùng khớp Gauss xuất hiện như là một lớp các phương pháp liên tục hấp dẫn để tích phân các PTVPCC giống như (1.15). Với ν = 1, phương pháp đã được biết đến như là phương pháp trung điểm, và với phương trình tổng quát (1.14) nó có dạng    h yn + y n+1 h yn+1 = yn + hf tn + , , x tn + − 1 . (1.17) 2 2 2 Áp dụng (1.17) cho (1.16) với cỡ bước tích phân hằng số h = 1/(m−δ), m ≥ 2, m nguyên, và 0 ≤ δ < 1, ta được 

yn + h a yn +y2 n+1 − ea π2 φ tn + h2 − 1

, tn + h2 − 1 ≤ 0, yn+1 = yn +yn+1 h h aπ yn + h a 2 −e 2η tn + 2 −1 , tn + 2 − 1 > 0,
trong đó η tn + h2 − 1 được cho bởi phép nội suy tuyến tính, đó là     η  = tn + 1 2 3 2 h −1 2 =η tn + h − (m − δ) h 2 − δ
yn−m + δ + 21 y
n−m+1 , 0 ≤ δ ≤ 12 − δ yn−m+1 + δ − 21 yn−m+2 , 12 < δ < 1

Tóm lại, công thức trung điểm cho (1.15) có dạng
a(n+ 1 )h 1 π π yn+1 = 1 + 2 ha yn − 2 he 2 sin 2 n+ 1 2


h−1 1 − 21 ha với n ≤ m − δ − 12 và  1 π a 1 1  (1+ 2 ha)yn − 2 e h(( 2 −δ1)yn−m +(δ+ 2 )yn−m+1 ) , 0 ≤ δ ≤ 1 , 2 1− ha yn+1 = 2 1 π a 3 1  (1+ 2 ha)yn − 2 e h(( 2 −δ)1yn−m+1 +(δ− 2 )yn−m+2 ) , 1− 2 ha 1 2 (1.18) (1.19) < δ < 1, với n > m − δ − 21 . Ta biết rằng, trong khi cấp chính xác rời rạc và đều của công thức trung điểm (ν = 1) cho các PTVPT đều bằng 2, thì công thức trùng khớp Gauss 2-nấc (ν = 2) lại có cấp chính xác rời rạc bằng 4 và cấp chính xác đều bằng 3. Để kiểm tra cấp chính xác rời rạc p của các phương pháp số giải PTVP với đầu ra liên tục thu được, xét nghiệm của (1.15) trong đoạn [0, 10] và nhớ rằng 20 nghiệm thuộc lớp C ∞ . Do đó không có điểm gián đoạn nào ảnh hưởng được đến cấp chính xác của phương pháp. Ta kí hiệu eh là sai số tuyệt đối lớn nhất tại các điểm nút với cỡ bước tích phân h = 1/(m − δ). Bằng cách chia đôi cỡ bước, giá trị tiệm cận của tỉ số rh = eh /e h được kì vọng bằng 2p . 2 Hình 1.8: Đồ thị logarit của các tỉ số rh , h = 1/(m − δ), m − δ = 2i ∗ k, như là một hàm số của i, cho các nghiệm số của phương trình (1.15) với a = 1 đạt được bằng phương pháp trùng khớp tại ν = 1 và ν = 2 các điểm Gauss. Các giá trị của m − δ được xác định bởi k = 1 (đường nét liền) và k = 5/3 (đường nét đứt). Với các giá trị nguyên (δ = 0) và không nguyên (δ > 0) của m − δ , tương ứng với h là hoặc không là “ước” của chậm τ = 1, công thức trung điểm duy trì được cấp chính xác rời rạc được kì vọng là 2. Ngược lại, công thức trùng khớp Gauss 2-nấc có cấp chính xác rời rạc là 4 hoặc 3 tuỳ theo m − δ là nguyên hay không nguyên, trong khi cấp chính xác đều là bằng 3. 1.3.2. Sự thất bại về tính ổn định Cũng giống như việc phân tích tính chính xác của các phương pháp số giải PTVP với đầu ra liên tục, ta xét lớp các phương trình thử tuyến tính hệ số hằng số sau  y 0 (t) = λy(t) − 45 λy(t − 1), t ≥ 0, y(t) = 1, t ≤ 0. (1.20) Phương trình này có nghiệm, được vẽ trên Hình 1.9 với một số giá trị của λ, là ổn định tiệm cận với mọi λ < 0. Công thức trung điểm (1.17), được mở rộng bởi phép nội suy tuyến tính, áp dụng cho phương trình (1.20) có dạng
1 4 yn+1 = 1 + 2 hλ yn − 5 hλ 1 − 12 hλ