Luận văn Về một phương pháp không cổ điển giải số phương trình vi phân bậc nhất và bậc hai

74 trang hapham 420

Download

Bạn đang xem 20 trang mẫu của tài liệu "Luận văn Về một phương pháp không cổ điển giải số phương trình vi phân bậc nhất và bậc hai", để 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:

luan_van_ve_mot_phuong_phap_khong_co_dien_giai_so_phuong_tri.pdf

Nội dung text: Luận văn Về một phương pháp không cổ điển giải số phương trình vi phân bậc nhất và bậc hai

BỘ GIÁO DỤC VÀ ĐÀO TẠO TRƯỜNG Về một phương pháp không cổ điển giải số phương trình vi phân bậc nhất và bậc hai 1
LỜI NÓI ĐẦU Phương trình vi phân là mô hình mô tả khá tốt các quá trình chuyển động trong tự nhiên và kĩ thuật. Để nghiên cứu phương trình vi phân, người ta thường tiếp cận theo hai hướng: nghiên cứu định tính và giải số. Mặc dù đã có lịch sử phát triển hàng trăm năm, do còn nhiều bài toán cần giải quyết, giải số phương trình vi phân thường vẫn thu hút sự quan tâm mạnh mẽ của các nhà toán học và các nhà nghiên cứu ứng dụng. Trong giải số phương trình vi phân, người ta thường cố gắng tìm ra những phương pháp hữu hiệu bảo đảm sự hội tụ, tính ổn định và tính chính xác cao. Để làm được điều này, người ta thường tổ hợp các phương pháp đa bước để nhận được các phương pháp mới có bậc hội tụ, tính ổn định và cấp chính xác cao hơn. Phương pháp không cổ điển giải số phương trình vi phân thường bậc nhất và bậc hai do M. V. Bulatov (và Berghe) đề xuất trong vòng năm năm trở lại đây nằm trong hướng này. Luận văn Về một phương pháp không cổ điển giải số phương trình vi phân bậc nhất và bậc hai có mục đích trình bày các phương pháp của Bulatov và Berghe theo các tài liệu [4] (2009) và [9]-[11] (2003-2008). Luận văn gồm ba Chương. Chương 1 trình bày một số khái niệm và phương pháp cơ bản giải số phương trình vi phân. Trong mục 1.2 của Chương, chúng tôi trình bày các phương pháp số cổ điển theo một quan điểm nhất quán là xuất phát từ Quy tắc cầu phương cơ bản. Chương 2 trình bày phương pháp không cổ điển (do Bulatov đề xuất vào những năm 2003-2008) giải số hệ phương trình vi phân bậc nhất, phi tuyến và tuyến tính, theo các tài liệu [9]-[11]. Chương 3 trình bày phương pháp không cổ điển giải số hệ phương trình vi phân bậc hai, tuyến tính và phi tuyến, theo bài báo của M. V. Bulatov và G. V. Berghe ([4], 2009). 2
Thông qua việc tính toán đạo hàm, phân tích các hàm nhiều biến vào chuỗi Taylor và các phép biến đổi chi tiết, chúng tôi cố gắng trình bày các kết quả của M. V. Bulatov và G. V. Berghe một cách rõ ràng và chi tiết nhất. Để minh họa và kiểm chứng lý thuyết, chúng tôi đã lập trình trên MATLAB và tính toán trên máy các ví dụ của M. V. Bulatov và G. V. Berghe. Luận văn được hoàn thành dưới sự hướng dẫn khoa học của PGS-TS Tạ Duy Phượng (Viện Toán học). Xin được tỏ lòng cám ơn chân thành nhất tới Thầy. Tác giả xin tỏ lòng cám ơn Ban Chủ nhiệm , các Thày Cô và các cán bộ khoa Toán- Cơ – Tin học, trường Đại học Khoa học Tự nhiên, Đại học Quốc gia Hà Nội đã nhiệt tình giảng dạy và giúp đỡ tác giả trong suốt quá trình học Cao học. Tác giả xin chân thành cám ơn Ban Lãnh đạo và các cán bộ, giáo viên Học viện Quân y đã tạo mọi điều kiện để tác giả hoàn thành tốt khóa học Cao học. Và cuối cùng, xin được cám ơn Gia đình, bạn bè đã thông cảm, sẻ chia, hy sinh và tạo mọi điều kiện cho tác giả trong suốt thời gian học Cao học và viết luận văn. Hà Nội, ngày 30 tháng 11 năm 2009 Tác giả Vũ Thị Thanh Bình 3
CHƯƠNG 1 Kiến thức chuẩn bị Trong Chương 1 chúng tôi nhắc lại những khái niệm cơ bản nhất của giải số phương trình vi phân nhằm thuận tiện cho trình bày ở các mục sau. 1.1. Bài toán Cauchy giải hệ phương trình vi phân Xét bài toán Cauchy tìm nghiệm của hệ phương trình x () t f ((),), x t t t  0,1 (1.1) thỏa mãn điều kiện ban đầu x(0) x0 , (1.2) trong đó f x,, t x t là các hàm vectơ n - chiều, hàm f xác định trên hình hộp vô hạn D : 0,1 R n . Ở đây ta hiểu nghiệm theo nghĩa cổ điển và địa phương, tức là nghiệm của (1.1)- (1.2) là một hàm khả vi x() t trên 0, , 1 sao cho x ( t ) f ( x ( t ), t ) trên 0, và x(0) x0 . Cùng với bài toán (1.1), ta cũng xét trường hợp hàm f(,) x t là tuyến tính, tức là f(,)()() x t B t x g t , trong đó B() t là ma trận cấp n n, còn g() t là vectơ n -chiều, tức là hệ tuyến tính x ( t ) B ( t ) x g ( t ), t  0,1. (1.3) Ta luôn giả thiết rằng các phần tử của ma trận B() t , của các vectơ f x, t , g() t là đủ trơn (có đạo hàm đến cấp cần thiết trong tính toán). Khi ấy theo định lí Picard- Lindelöf, hệ (1.1)-(1.2) có duy nhất nghiệm x() t trên toàn đoạn 0,1 (nghiệm có thể kéo dài được trên toàn bộ khoảng xác định, hay tồn tại nghiệm toàn cục, xem [8], trang 467). Lưu ý này là quan trọng trong giải số hệ phương trình (1.1)-(1.2). 4
1.2. Giải số bài toán Cauchy Để chứng minh định lý về sự tồn tại và duy nhất nghiệm của hệ phương trình vi phân (1.1)-(1.2), ta có thể xây dựng dãy nghiệm xấp xỉ hội tụ tới nghiệm của bài toán (1.1)-(1.2) trên khoảng tồn tại nghiệm. Có hai phương pháp xây dựng dãy nghiệm xấp xỉ: phương pháp giải tích và phương pháp số kết quả được cho dưới dạng bảng, như phương pháp Euler, phương pháp Runge-Kutta, phương pháp đa bước, Dưới đây trình bày cách xây dựng các công thức Euler, Runge-Kutta, xuất phát từ qui tắc cầu phương cơ bản (xem, thí dụ, [2]). 1.2.1. Quy tắc cầu phương cơ bản và giải số phương trình vi phân Quy tắc cầu phương cơ bản (basic quadrature rules) có thể được coi là phương pháp quan trọng để tính tích phân. Vì giải phương trình vi phân thường (1.1) với điều kiện ban đầu (1.2) tương đương với việc giải phương trình tích phân t x( t ) x f ( x ( s ), s ) ds (1.4) 0 t0 nên ta cũng có thể sử dụng quy tắc cầu phương cơ bản trong việc giải số phương trình vi phân. Trong mục này ta sẽ chỉ ra rằng, nhiều công thức sai phân cổ điển giải số phương trình vi phân có thể suy ra từ quy tắc cầu phương cơ bản. Trước tiên ta nhắc lại quy tắc cầu phương cơ bản (xem, thí dụ, [1]). b Nội dung cơ bản của quy tắc cầu phương là: để tính tích phân f() t dt ta thay f() t a bởi một đa thức nội suy (interpolating polynomial). Tích phân của hàm f() t được xấp xỉ bởi tích phân của hàm đa thức (tính được chính xác). Giả sử ta có s điểm nội suy khác nhau c1, c 2 , , cs trong khoảng a, b . Đa thức nội suy Lagrange bậc nhỏ hơn s có dạng (xem [1]): s ()()()t  f cj L j t , j 1 5
s ()t c b s trong đó L() t i . Khi ấy f()() t dt  f c . j   j j i 1, i j ()cj c i a j 1 b Các trọng số  được tính theo công thức  L(). t dt j j j a Nếu s 1 thì đa thức nội suy ()()t f c1 và ta có: b f( t ) dt ( b a ) f ( c ). 1 a Ta nói độ chính xác (precision) của quy tắc cầu phương là p nếu quy tắc này chính xác cho mọi đa thức bậc nhỏ hơn p , tức là với mọi đa thức Pk () t bậc nhỏ hơn p ta có: b s P( t ) dt  f ( c ). k j j a j 1 Nếu b a 0( h ) thì sai số trong quy tắc cầu phương của độ chính xác p là 0(h p 1 ). Ta xét một số trường hợp đặc biệt. Nếu chọn s 1 và c1 a thì ta có công thức xấp xỉ tích phân bởi diện tích hình chữ nhật ABCD (Hình 1.1): b f( t ) dt ( b a ) f ( a ). (1.5) a Nếu x() t là nghiệm của phương trình vi phân (1.1) - (1.2) (nghiệm của phương trình tích phân (1.4)) thì: t h xth( ) xt ( ) fxssds ( ( ), ) (1.6) t Kết hợp với công thức (1.5) ta đi đến công thức: x( t h ) x ( t ) h . f ( x ( t ), t ) (1.7) Gọi h là độ dài bước (stepsize) của biến độc lập t ( h có thể dương hoặc âm, khi h dương thì nghiệm được xây dựng về bên phải của điểm t và ngược lại, khi h âm 0 6
thì nghiệm được xây dựng về bên trái của t0 ). Dưới đây ta coi h 0 , trường hợp h 0 có thể được xét tương tự. f Từ công thức (1.7) ta có x x h.(,) f x t ; 1 0 0 0 D C x2 x 1 h.(,) f x 1 t 1 ; ; a b O A B x xn 1 x n h.(,) f x n t n . Đây chính là công thức Euler tiến quen thuộc. Hình 1.1 Nếu chọn s 1 và c1 b thì ta có công thức xấp xỉ tích phân bởi diện tích hình chữ nhật ABEF (Hình 1.2): f F E b f( t ) dt ( b a ) f ( b ). a Từ đây ta có: a b xth( ) xt ( ) hfxthth . ( ( ), ) O A B x Suy ra công thức Euler lùi: xn 1 x n h.(,) f x n 1 t n 1 . Hình 1.2 Hai phương pháp Euler tiến và Euler lùi là những phương pháp Runge-Kutta bậc nhất (có độ xấp xỉ bậc nhất). f a b Nếu chọn s 1 và c thì ta có 1 2 N M công thức xấp xỉ tích phân bởi diện tích hình chữ nhật ABMN (Hình 1.3): a b t h x h h O A a b B fxssds((),) hfxt .(( ), t ) 2 2 2 t Từ đây ta có: Hình 1.3 h h x( t h ) x ( t ) h . f ( x ( t ), t ) . 2 2 7
Từ công thức trên ta có h h x x h[ f ( x ( t ), t ] . n 1 n n2 n 2 Đây chính là phương pháp trung điểm (midpoind method). t b t a Nếu chọn s 2 và c a, c b thì L() t và L(). t 1 2 1 ()a b 2 ()b a Suy ra b b b t b1 ( t b )2 b a  L() t dt dt 1 1 a a (a b ) ( a b ) 2a 2 và a b b t a1 ( t a )2 b a  L(). t dt dt 2 2 a a (b a ) ( b a ) 2b 2 Chứng tỏ b b a f( t ) dt [ f ( a ) f ( b )]. a 2 t h Như vậy nếu xấp xỉ tích phân f( x ( s ), s ) ds bởi công thức trên (bởi diện tích hình t thang ABED, Hình 1.4) thì ta được: f E t h h fxssds((),) [(( fxthth ), ) fxtt ((),)]. t 2 D Từ đây ta có công thức hình thang: h b x x [ f ( x , t ) f ( x , t )]. a n 1 n2 n n n 1 n 1 O A B x Phương pháp điểm giữa và phương pháp hình thang là hai phương pháp ẩn, Hình 1.4 chúng có độ chính xác p 2 . a b Nếu chọn s 3 và c a,, c c b thì, đặt h b a , ta có: 1 22 3 8
a b (t )( t b ) 2 a b L( t ) 2 ( t )( t b ), 1 a b 2 (a )( a b ) h 2 2 (t a )( t b ) 4 L( t ) ( t a )( t b ), 2 a b a b 2 ( a )( b ) h 2 2 a b (t a )( t ) 2 a b L( t ) 2 ( t a )( t ). 3 a b 2 (b a )( b ) h 2 2 Suy ra b2 ba b 2 b a b  Ltdt( ) ( t )( tbdt ) ( tb )( tbdt ) 1 1 2 2 ah a2 h a 2 b 2b a b 2 ( t b )3 a b ( t b ) 2 [()()()]t b2 t b dt [ () ] 2 2 2 3 2 2 ha h a 2 (b a )3 h . h2 12 6 và b 4 b 4 b  Ltdt() ()() tatbdt ()( tata ()) abdt 2 2 2 2 ah a h a b 4 (t a )3 ( t a ) 2 4 h [ (a b )] . 2 3 2 6 h a Do tính chất đối xứng (hoặc tính trực tiếp), ta có h   . 3 1 6 Từ các tính toán trên ta đi đến công thức Simpson: b h a b f() t dt [()4( f a f ) f ()] b . a 6 2 Suy ra công thức xấp xỉ nghiệm của phương trình vi phân 9
h h h xthxt( ) () [((),)4(( fxtt fxt ), t ) fxthth (( ), )] 6 2 2 và công thức sai phân h h h xx [(,)4(( fxt fxt ), t ) fxt ( , )] . n 1 n6 n n n 2 n 2 n 1 n 1 Đây là công thức ẩn của phương pháp Runge-Kutta kinh điển cấp bốn (classical fourth-order Runge-Kutta method). 1.2.2. Phương pháp Runge-Kutta 1.2.2.1. Dẫn tới phương pháp Runge - Kutta Vì phương pháp ẩn đòi hỏi tại mỗi bước phải giải một phương trình phi tuyến, điều này không đơn giản, nên ta cố gắng xây dựng các công thức Runge-Kutta hiển từ công thức hình thang ẩn, công thức điểm giữa ẩn và công thức Runge-Kutta kinh điển cấp bốn ẩn tương ứng như sau. Trong công thức hình thang ẩn: h x x [ f ( x , t ) f ( x , t )], n 1 n2 n n n 1 n 1 ta thay giá trị xn 1 ở vế phải bằng công thức Euler tiến: xˆn 1 x n hf(,) x n t n . Khi ấy ta được công thức: h x x [ f ( x , t ) f ( xˆ , t )] n 1 n2 n n n 1 n 1 Công thức này được gọi là phương pháp hình thang hiển (explicit trapezoidal method). h Bằng cách sử dụng xấp xỉ bậc nhất của x() t theo phương pháp Euler tiến: n 2 h xˆ 1 xn f(,) x n t n n 2 2 và thay vào công thức của phương pháp trung điểm ẩn h h x x hf((,) x t t n 1 n n2 n 2 10
ta nhận được phương pháp trung điểm hiển (explicit midpoint method): h xn 1 x n hf(,) xˆ 1 t n n 2 2 Từ phương pháp Runge-Kutta ẩn cấp bốn kinh điển h h h xxn 1 n [(,)4(( fxt n n fxt n ), t n ) fxt ( n 1 , n 1 )] , 6 2 2 ta có công thức Runge-Kutta hiển bậc bốn kinh điển sau: h x x ( k 2 k 2 k k ), n 0,1,2, n 1 n 6 1 2 3 4 trong đó: k1 f( xn , t n ); hk1 h k2 f( xn , t n ); 2 2 hk h k f( x 2 , t ); 3 n2 n 2 k4 f( xn hk 3 , t n ). 1.2.2.2. Phương pháp Runge-Kutta tổng quát Nội dung cơ bản của phương pháp Runge-Kutta tổng quát như sau. Chia đoạn 0,1 thành một lưới đều 1 t ih , i 0,1,2, , N , h , i N và kí hiệu xi là giá trị xấp xỉ x ti , Bi B() t i , gi g t i . Phương pháp Runge-Kutta cho bài toán (1.1)-(1.2) có dạng (xem, [2], [4]-[7]) s xi 1 x i h b i X i , (1.8) i 1 trong đó X i là vectơ n -chiều và là nghiệm của hệ phương trình phi tuyến s X f x h a X , t c h . (1.9) i i  i j j i i j 1 Các tham số aij , ci , bi xác định bậc xấp xỉ của phương pháp, còn s được gọi là số 11
nấc. Nếu aij 0 với mọi j i thì ta có phương pháp Runge-Kutta hiển. Khi ấy tính toán khá đơn giản ( X i được tính theo công thức truy hồi). Nếu aij 0 với j i nào đó thì ta có phương pháp Runge-Kutta ẩn. Khi ấy tại mỗi bước ta phải giải một hệ ns phương trình phi tuyến (tuyến tính nếu f(,)()() x t B t x g t ) để tìm s vectơ X i (mỗi vectơ X i có n tọa độ). Thường phương pháp Runge-Kutta được viết dưới dạng bảng Butcher (Butcher table) c A bT Hai phương pháp Runge-Kutta quan trọng thường hay được sử dụng là phương pháp Runge-Kutta bậc hai và phương pháp Runge-Kutta bậc bốn. 1.2.2.3. Công thức lặp của phương pháp Runge-Kutta bậc hai Giả thiết rằng ta đã biết giá trị của x tại tn là xn . Phương pháp Runge-Kutta hiển hai nấc cấp hai sử dụng điểm (,)x t để xấp xỉ giá trị của x tại điểm tiếp theo bằng n n công thức xn 1 x n h( b 1 k 1 b 2 k 2 ), (1.10) trong đó k1 fxy(n , n ); k 2 fx ( n chy 2 , n hak 21 1 ). Khái niệm s -nấc ( s -stage) thể hiện rằng số lần tính các giá trị của hàm f (tại các điểm khác nhau trong công thức Runge-Kutta) là s . Để tìm các phương pháp Runge-Kutta bậc hai, ta làm như sau (xem [2]). Khai triển Taylor hàm f(,) x t theo phương trình (1.1) và theo công thức (1.10) rồi so sánh, ta đi đến kết luận: Các hệ số trong phương pháp Runge-Kutta cấp hai phải thoả mãn hệ phương trình 1 1 b b 1, c b , a b . 1 2 2 22 21 2 2 12
Đây là một hệ ba phương trình (phi tuyến) bốn ẩn. Ta có thể chọn một hệ số, thí dụ, b2 0 tự do. Khi ấy các hệ số còn lại biểu diễn qua b2 0 bởi các công thức: 1 1 b1 1 b 2 , c2 , a21 . 2b2 2b2 1 1 Chọn b , thì b và a c 1. Khi ấy ta có một phương pháp Runge-Kutta 2 2 1 2 21 2 cấp hai cho phép tính xn 1 dựa trên công thức: 1 1 x x hfxt(,)((,)) hfx hfxtt h . n 1 n2 n n 2 n n n n Công thức này được gọi là phương pháp Runge-Kutta đơn giản (Simple Runge- Kutta Method) hoặc phương pháp tiếp tuyến cải tiến (Impoved Tangent Method), vì nó trùng với phương pháp Euler cải tiến. 1 1 Nếu chọn b 1 thì a , b 0 và c . Khi ấy ta có công thức 2 21 2 1 2 2 h h x xhbkbk ( ) x . hfx . ( fxtt ( , ), ) . n 1 n 1 1 2 2 n n2 n n n 2 Phương pháp tính theo công thức trên được gọi là phương pháp Euler-Cauchy. 1.2.3. Phương pháp cổ điển đa bước Phương pháp cổ điển k -bước cho bài toán (1.1) có dạng (xem, [3], [4]-[7]) k k  jx i 1 j h   j f(,) x i 1 j t i 1 j . (1.11) j 0 j 0 Phương pháp một tựa tương ứng với nó là k k k  jx i 1 j hf(,)   j x i 1 j   j t i 1 j . (1.12) j 0 j 0 j 0 Đối với phương pháp này ta giả thiết rằng các giá trị xuất phát x1, x 2 , , xk 1 đã được tính tương đối chính xác. Nếu 0 0 và  0 0 thì phương pháp là phương pháp hiển. Nếu 0 0 và  0 0 thì ta có phương pháp ẩn. 13
1.3. Mô hình thử và ổn định của phương pháp số 1.3.1. Mô hình thử Để phân tích hiệu quả của các phương pháp, ta thường thử chúng trên mô hình G. Dahlquist (gọi là phương trình thử hay mô hình thử) x  x, x 0 x0 R , t  0,1, (1.13) trong đó  là một hằng số (thực hoặc phức). Nghiệm của phương trình này là t x(t) e x0 . Ta thường viết  R iI , trong đó R và I tương ứng là các phần thực và phần ảo của  . Tương ứng với phương trình (1.13), xét phương trình sai phân tuyến tính bậc nhất xn 1 xn , n 0,1,2, , trong đó x0 cho trước và  nói chung là một số phức. Nghiệm của phương trình n này là xn  x0 . Ta thấy rằng nghiệm này bị chặn khi và chỉ khi  1. Giả sử bước h 0 cố định. Khi ấy giá trị của nghiệm chính xác tại các điểm tn nh tn nh n h sẽ là xn e x0 e x0  x0 , trong đó  e . Nếu nghiệm chính xác bị chặn thì  eh 1. Điều này chỉ có thể xảy ra nếu Re h R h 0 . Điều này có nghĩa là, trên mặt phẳng với trục hoành Re(h) và trục tung Im(h) , miền ổn định của nghiệm chính xác phải là nửa mặt phẳng mở bên trái. Phương pháp một bước được gọi là ổn định tuyệt đối nếu  1 và ổn định tương đối nếu  e h . Nếu  là thuần ảo và  1thì ổn định tuyệt đối được gọi là ổn định tuần hoàn (P- ổn định). Khi miền ổn định của phương trình sai phân đồng nhất với miền ổn định của phương trình vi phân, lược đồ sai phân hữu hạn được gọi là ổn định - A. 14
Phương trình thử thường được sử dụng như một mô hình để dự đoán tính ổn định của phương pháp số giải hệ dạng tổng quát (1.1)-(1.2). Để thuận tiện, ta cũng có thể đưa ra các khái niệm ổn định tương tự như sau. Kí hiệu z  h , mọi phương pháp Runge-Kutta (1.8)-(1.9) đều có thể viết dưới dạng xi 1 R() z x i , trong đó R() z được gọi là hàm ổn định. Định nghĩa 3.1 Tập tất cả các điểm của mặt phẳng phức M mà R( z ) 1 được gọi là miền ổn định của phương pháp (1.8)-(1.9). Nếu tập hợp đó chứa toàn bộ nửa mặt phẳng trái thì phương pháp được gọi là ổn định-A, còn nếu ngoài ra limR ( z ) 0 thì phương pháp z được gọi là ổn định-L(hay còn gọi là ổn định tiệm cận). 1.3.2. Sự ổn định của phương pháp Euler Phương pháp Euler áp dụng cho phương trình thử (3.1) có dạng xn 1 xn hf (xn ,tn ) xn hxn 1 h xn . Nghiệm của phương trình sai phân tương ứng là n n xn 1 h xn  x0 , trong đó  1 h . Phương pháp số là ổn định nếu  1. Xét các trường hợp sau 1)  là số thực. Khi ấy  1 h 1, hay 2 h 0. 2)  là thuần ảo ( i , trong đó  là số thực khác 0). Khi ấy  1 ih 1  2 h 2 1. Chứng tỏ phương pháp là không ổn định nếu  là thuần ảo. 3)  là số phức ( R iI ). Khi ấy 2 2  1 h 1 R h iI h 1 R h I h 1, 15
nghĩa là h nằm trong hình tròn đơn vị tâm là 1;0 (Hình1.5). Hình tròn này tiếp xúc với trục ảo. Im h - -1 Re  O h Hình 1.5 Như vậy, chỉ có một phần rất nhỏ (hình tròn bán kính bằng 1) của nửa mặt phẳng trái là miền ổn định của phương pháp Euler. Với mọi giá trị khác của h trong nửa mặt phẳng trái và bên ngoài hình tròn này, nghiệm số sẽ bị phóng đại (blow-up) khi nghiệm chính xác triệt tiêu (decays). Phương pháp số này được gọi là ổn định có điều kiện. Để nhận được nghiệm số ổn định, bước h phải được chọn sao cho h nằm trong 2 hình tròn. Nếu  là số thực âm thì từ điều kiện 2 h 0 suy ra h 0 .  Nếu  là thực và nghiệm số không ổn định thì 1 h 1, nghĩa là 1 h là một số n âm và có trị tuyệt đối lớn hơn 1. Vì xn 1 h x0 nên nghiệm số sẽ đổi dấu qua mỗi bước. Sự thay đổi của nghiệm số mô tả khá rõ tính không ổn định. Tương tự, ta có thể xét tính ổn định của các phương pháp Euler cải tiến. Ta đi đến kết luận sau. Phương pháp Euler ẩn là ổn định-L và hội tụ cấp một; Phương pháp hình thang là ổn định-A và hội tụ cấp hai, còn phương pháp Euler hiển không phải là ổn định-A và hội tụ cấp 1. z 1 1 Hàm ổn định của các phương pháp này tương ứng là , 2 và 1 z . z 1 z 1 2 16
1.3.3. Sự ổn định của phương pháp Runge-Kutta 1.3.3.1 Sự ổn định của phương pháp Runge-Kutta bậc hai Xét phương pháp Runge-Kutta bậc hai cho phương trình thử (1.13). Ta có k1 hf (xn ,tn ) hxn ; k2 hf (xn k1,tn h) h xn k1 h xn hxn h 1 h xn và 1 1 2h 2 xn 1 xn k1 k 2 xn h h 1 h xn 1 h xn . 2 2 2 2h 2 Để phương pháp ổn định thì  1, trong đó  1 h . 2 2h 2 Trường hợp 1.  là số thực. Khi ấy 1 h 1 hay 2 h 0. 2 Trường hơp 2.  i thuần ảo,  0 . 1 Khi ấy  1  4h 4 1. Phương pháp không ổn định. 4 Trường hợp 3.  R iI là số phức. Khi ấy 2 h2 2 2 1 h 1  h R I h 2 i    R I R I 2 2 2h 2 là số phức. Đặt 1 h ei và tìm nghiệm phức h của phương trình bậc hai 2 theo các giá trị của  . Nhận xét rằng  1 với mọi giá trị của  . Miền ổn định được chỉ ra trên Hình 1.6. 1.3.3.2. Sự ổn định của phương pháp Runge-Kutta bậc bốn Xét phương pháp Runge-Kutta bậc bốn cho phương trình thử (1.13). Ta có k1 hf (xn ,tn ) hxn ; k1 h k2 h xn h xn hxn h 1 xn ; 2 2 17
k 1 h h 2 h 2 2 k3 h xn h xn h 1 xn h 1 xn ; 2 2 2 2 4 h 2 h2 2h 2 3h3 k h x k h x h 1 x h 1 h x 4 n 3 n n n 2 4 2 4 Và 1 2h 2 3h 3 4h 4 xn 1 xn k1 2k2 2k3 k 4 1 h xn . 6 2 3! 4! 2 h2 3h3 4 h4 Để phương pháp ổn định thì  1, trong đó  1 h . 2 3! 4! Trường hợp 1.  là số thực. Khi ấy 2.785 h 0. Trường hơp 2.  i thuần ảo,  0 . Khi ấy 0 h 2 2 . 2 h 2 3h3 4h 4 Trường hợp 3.   i là số phức. Đặt 1 h ei và tìm R I 2 3! 4! nghiệm phức h của phương trình bậc bốn theo các giá trị của  . Nhận xét rằng  1 với mọi giá trị của  . Miền ổn định được chỉ ra trên Hình 1.6. Hình 1.6 18
1.3.4. Sự ổn định của phương pháp đa bước Áp dụng các phương pháp (1.11) và (1.12) cho bài toán (1.13) ta được k  j z  j x i 1 j 0 . (1.14) j 0 Phương trình đặc trưng của phương trình sai phân tuyến tính trên có dạng k k j  j z  j  0. (1.15) j 0 Định nghĩa 3.2 Tập tất cả các điểm của mặt phẳng phức M mà  (z ) 1 với mọi nghiệm của (1.15) và  (z ) 1 đối với các nghiệm bội được gọi là miền ổn định của phương pháp (1.14). Nếu tập hợp đó chứa toàn bộ nửa mặt phẳng trái thì phương pháp được gọi là ổn định-A. Nhận xét Với mọi phương pháp cổ điển miền ổn định chứa gốc tọa độ của mặt phẳng phức. Tương ứng với bậc của phương pháp ta đã biết những điều sau (xem [11]): 1) Bậc của các phương pháp Runge-Kutta s -nấc cho phương trình (1.13)-(1.15) không vượt quá 2s (chắn Butcher). 2) p là cấp chính xác, k là số bước của phương pháp (1.11).Nếu phương pháp ổn định thì p không vượt quá k 2 khi k chẵn và k 1 khi k lẻ (chắn Dahlquist thứ nhất). 3) Phương pháp (1.11) ổn định – A không thể có cấp chính xác vượt quá 2 (chắn Dahlquist thứ hai). Trong Chương sau ta sẽ trình bày phương pháp do Bulatov đề nghị cải tiến được những hạn chế nêu trên. 1.3.5. Sự ổn định của phương pháp sai phân hữu hạn Xét phương trình vi phân tuyến tính bậc hai x kx 0, trong đó k là hằng số và đủ lớn so với 1. 19
Phương trình có nghiệm là x(t) A Be kx , trong đó A và B là các hằng số bất kì được xác định bởi điều kiện ban đầu hoặc điều kiện biên tương ứng. Nếu k 0 và x thì nghiệm bị chặn. Đại lượng e kx là nghịch biến khi k 0 và đồng biến khi k 0. Xấp xỉ sai phân trung tâm cho hệ là 1 k x 2x x x x 0 . h 2 i 1 i i 1 2h i 1 i 1 i 2 kh Phương trình sai phân này có nghiệm là xi A B . 2 kh i 2 kh 2 kh Nếu k 0 đại lượng 1 nên là nghịch biến. Khi k 0 thì 2 kh 2 kh 2 kh 2 1 nếu h . Đây chính là điều kiện ổn định cho hệ sai phân. 2 kh k 20
CHƯƠNG 2 Về một phương pháp không cổ điển giải số hệ phương trình vi phân cấp một Chương này trình bày một phương pháp mới do Bulatov đề xuất giải số bài toán Cauchy cho hệ phương trình vi phân cấp một (xem [9]-[11]) tốt hơn các phương pháp cổ điển. Phương pháp mới là một họ phương pháp một bước, bậc hai, trong đó có phương pháp là L-ổn định. Nội dung của Chương gồm hai mục. Trong 2.1 chúng tôi trình bày phương pháp không cổ điển do Bulatov đề xuất giải số hệ phương trình vi phân phi tuyến cấp một. Phương pháp không cổ điển do Bulatov đề xuất giải số hệ phương trình vi phân tuyến tính cấp một được trình bày trong 2.2. Để làm sáng tỏ phương pháp, chúng tôi đã thực hiện các tính toán chi tiết (phân tích các hàm nhiều biến dưới dạng chuỗi Taylor, ) mà trong [9]-[11] trình bày không tường minh. 2.1. Phương pháp không cổ điển giải số hệ phương trình vi phân phi tuyến cấp một 2.1.1. Phương pháp tổng quát Xét bài toán Cauchy tìm nghiệm của hệ phương trình x () t f ((),), x t t t  0,1 (2.1) thỏa mãn điều kiện ban đầu x(0) x0 , (2.2) trong đó f x,, t x t là các hàm vectơ n - chiều, hàm f xác định trên hình hộp chữ nhật vô tận D : R n 0,1. Để giải bài toán (2.1)-(2.2), ta bắt đầu đi từ phương pháp- và phương pháp một tựa của nó (hai phương pháp này đều có cấp chính xác bằng một): 21
xi 1 xhcfxt i  1( i 1 , i 1 ) (1 cfxt 1 ) ( i , i  (2.3) và xi 1 xhfcx i ( 2 i 1 (1 cxtch 2 ) i , i 2 ) (1.4) với c1, c 2 là các hằng số tùy ý, 0 c1 , c 2 1. Mỗi công thức truy hồi (2.3) (hoặc (2.4)) cho một dãy các giá trị xấp xỉ nghiệm của phương trình vi phân (2.1)-(2.2). Dưới đây ta cố gắng kết hợp hai phương pháp (2.3) và (2.4) để được một phương pháp số mới giải hệ phương trình vi phân (2.1)- (2.2). Khai triển Taylor tại điểm x xi ta được: 2 fxt i 1,,, i 1 fxt i i 1 Jxt i i 1 x i 1 xOx i i 1 x i và f c2 xi 1 1 c 2 x i , t i c 2 h 2 fxtchJxtchcx i, i 2 i , i 2 2 i 1 1 cx 2 i x i Ox i 1 x i fxtch,,, cJxtchx x Ox x 2 i i2 2 i i 2 i 1 i i 1 i f x, t trong đó J x, t . x Các phương pháp (2.3)-(2.4) được tuyến tính hóa như sau xi 1 xhcfxt i  1(,)(1 i i 1 cfxt 1 )(,) i i cJxt 1 (,)( i i 1 x i 1 x i ) , (2.3’) xi 1 x i hfxt ( i , i ch 2 ) cJxt 2 ( i , i chx 2 )( i 1 x i )] . (2.4’) Viết lại hai công thức trên thành một hệ phương trình đại số tuyến tính có ẩn là xi 1 (kí hiệu E là ma trận đơn vị cấp n ): EhcJxt 1(,)i i 1 EhcJxt 1 (,) i i 1 cfxt 1 (,)(1)(,) i i 1 cfxt 1 i i xi 1 x i h . (2.5) EhcJxt 2(,)(,)(,)i i ch 2 EhcJxt 2 i i ch 2 cfxt 2 i i ch 2 Nhận thấy rằng hệ (2.5) nói chung không có nghiệm theo nghĩa cổ điển vì số phương trình nhiều hơn số ẩn, tức là hệ (2.3) và (2.4) nói chung không có nghiệm trùng nhau. Để giải hệ phương trình đại số (2.5) ta nhân hai vế của hệ này với ma 22
trận E hcJxt1(,)((,)i i 1 cE 3 hcJxt 2 i i ch 2 cấp n 2 n , trong đó c3 là hằng số tùy ý. Khi ấy (2.5) trở thành: 2 2 EhcJxt ,, cEhcJxt ch x 1i i 1 3 2 i i 2 i 1 2 2 EhcJxt ,, cEhcJxt ch x 1i i 1 3 2 i i 2 i (2.6) hEhcJxt 1 i, i 1 cfxt 1 i , i 1 1 cfxt 1 i , i hc3 E hcJxt 2 i,,. i ch 2 fxt i i ch 2 Theo [11]lược đồ sai phân (2.6) là ổn định với mọi bộ hệ số c1,, c 2 c 3 và có bậc hội tụ tối thiểu là bậc một. Ta có ba tham số tự do, vì vậy có thể chọn được bộ ba số c1,, c 2 c 3 sao cho họ lược đồ sai phân (2.6) hội tụ cấp hai. Khai triển Taylor theo t tại điểm xi, t i ta được: 2 Jxt i,,, i 1 Jxt i i JxthOh t i i 2 Jxt i,,, i ch2 Jxt i i JxtchOh t i i 2 2 fxt i,,, i 1 fxt i i fxthOh t i i 2 fxt i,,, i ch2 fxt i i fxtchOh t i i 2 h2 x x t x hx x O h3 i 1 i 1 i i2 i Thay vào (2.6) ta được (1) cEh3 2(,)2 cJxt 1i i chJxt 1 t i ,2 i ccJxt 2 3 (,)2 i i cchJxt 2 3 t i , i 2 2 2 h hcJxt2 2,,,, hJxt ccJxt 2 chJxt xhx x 1 iitii 2 3 ii 2 tii ii i  2 (1) cEh 2(,)2 cJxt chJxt , 2 ccJxt (,)2 cchJ x, t 3 1i i 1 t i i 2 3 i i 2 3 t i i 2 2 hcJxt2 2 ,,,, hJxt ccJxt 2 chJxt x 1 ii tii 3 2 ii 2 tii  i hEhcJxt 1(ii , ) hcJxt 1 tii , cfxt 1 ii , chfxt 1 tii , 1 cfxt 1 ii , 2 2 3 hcEhcJxt3 2(,),,,i i hcJxt 2 t i i fxt i i hcfxt 2 t i i Oh hay 23
h2 (1 cE ) 2 hcccJxt , h2 2 cccJxt , cccJxt 2 2 2 , hx x 3 1 2 3 i i 1 2 3 t i i 1 2 3 i i  i i 2 2 2 3 hEhcJxt 1 ii,,,,,, fxt ii hcfxt 1 tii hcEhccJxt 3 2 3 ii fxt ii hcfxt 2 tii Oh hay h2 (1 cE ) 2 hcccJxt , h2 2 cccJxt , c 2 ccJxt 2 2 , hx x 3 1 2 3 i i 1 2 3 t i i 1 2 3 i i  i i 2 2 2 2 3 hfxt ii,,,,,,, chfxt1 tii chJxt 1 ii fxt ii chJxt 1 iitii fxt hcfxt 3 ii 2 2 2 3 3 cchf2 3t xt i,,,, i cchJxt 2 3 i i fxt i i cchJxt 2 3 i i ft x i,. t i O h Đẳng thức trên đúng với mọi h . So sánh hệ số hai vế ta được: Hệ số của h1 : (1 c3 )xi f (xi ,ti ) c3 f (xi ,ti ) . Do x ( t ) f ( x ( t ), t ) (phương trình (1.1)) nên hệ thức trên đúng với mọi c3 . Hệ số của h2 : x (1 c )i 2( c c c ) J ( x , t ) x 32 1 2 3 i i i cfxt1tii (,) cJxtfxt 1 (,)(,) ii ii ccfxt 2 3 tii (,) ccfxtJxt 2 3 (,)(,). iiii Từ phương trình (1.1) ta có: x ( t ) f ( x ( t ), t ) và f  f x x Jxtfxt ,,, fxt . tx t  t t Do đó hệ số của h2 viết lại thành: 1 c 3 (Jxtfxt , , fxt (,))( cccJxtfxt )((,)(,) fxt ,) 2 ii iitii1 2 3 iiiitii hay 1 c3 c1 cc 2 3 Jxtfxt i, i i , i fxt t ( i , i ) 0 . 2 Như vậy, hệ số của h2 bằng 0 khi và chỉ khi ta có 1 c3 2 c 1 2 c 2 c 3 . (2.7) 24
Như vậy nếu chọn c1,, c 2 c 3 thỏa mãn (2.7) thì ta được công thức (2.6) có cấp hai, bởi vì lúc này công thức (2.6) có sai số địa phương bậc h3 . Ta đi đến định lý sau. Định lý 1.1 2 2 Nếu f x t ,t , J x t ,t C 0,1 thì ta có đánh giá xi x t i O() h , trong đó xi tìm được theo công thức (2.6). Định lý này được chứng minh nhờ nhận xét là sai số địa phương có bậc ba, và ma trận chuyển sang bước tiếp theo có chuẩn là 1 hK với 0 K . 2.1.2. Phương trình thử Áp dụng công thức (2.6) cho phương trình thử x  x, x 0 x0 R , t  0,1 . (2.8) Ta có f x,, t  x do đó J x,,, t   x  t và f(,),(,) xi t i c2 h  x i f x i 1 t i 1  x i 1 . Khi đó (1.6) được viết lại thành: Ehc  2 cEhc  2 x Ehc  2 cEhc  2 x 1 3 2 i 1 1 3 2 i hEhc( 1 )( cx 1 i (1 c 1 )  x i ) hcE 3 (  hc 2 )  x i hay (1 cE ) (2 c 2 cch ) ( ccch2 2 )()]  2 x 3 1 2 3 1 2 3i 1 (2.9) 2 2 2 (1)(1 cE3 c 3123 22) c cch ( cccccch 123123 )()].  xi Đặt z h và nhận xét rằng, với phương trình thử, ma trận E 1 là một số, ta được công thức xi 1 R( z ). x i , trong đó 2 2 2 (1 c3 ) (1 2 cc 13 2 cczccccccz 23 ) ( 113223 ) R z 2 2 2 . (1 c3 ) 2( c 1 c 2 c 3 ) z ( c 1 c 2 c 3 ) z Để (2.6) là ổn định-L thì limR z 0 , suy ra trong R z bậc của tử thức phải nhỏ z hơn bậc của mẫu thức, tức là ta phải chọn c1,, c 2 c 3 thỏa mãn 25
2 2 c1 c 2 c 3 0, 2 2 c1 c 1 c 3 c 2 c 2 c 3 0, (2.10) c 1. 3 Chú ý rằng khi c3 1 thì điều kiện (2.7) cho c1 c3c2 0 . Khi ấy R z 1 2.1.3. Trường hợp đặc biệt Nếu chọn c1 1, c 2 0, c 3 1 thì công thức (2.6) trở thành 2 E E hJ xi, t i 1 x i 1 (2.11) xi hfxt i,,,,. i EhJxt i i 1 EhJxt i i 1 x i hfxt i i 1 1 Hàm ổn định của (2.11) cho phương trình thử là R() z có limR ( z ) 0 . z2 z 1 z 2 Như vậy (2.11) là công thức có tính chất ổn định-L. Trường hợp đặc biệt này đã được xét độc lập trong [9]. Công thức (2.11) trùng với một phương pháp Runge-Kutta, đó là phương pháp Lobatto III C với bảng Butcher 0 ½ -1/2 1 ½ ½ ½ ½ Tuy nhiên, trong trường hợp tuyến tính, khi ma trận B() t phụ thuộc vào t , thì các phương pháp (2.11) và phương pháp Lobatto cho các kết quả khác nhau. Nhận xét Khi thực hiện phương pháp Lobatto III C cho bài toán (1.1)-(1.2), tại mỗi bước ta cần giải hệ phương trình phi tuyến cấp 2n 2 n . Những hệ này thường được giải theo phương pháp Newton cải tiến, trong đó gần đúng ban đầu được chọn bởi giá trị đã tính được ở bước trước. Để thực hiện được phương pháp, ta cần phải đưa ra tiêu chuẩn chọn số bước lặp. Các phương pháp giải hệ phương trình tuyến tính cấp n3 n n đòi hỏi khoảng phép toán số học và phép nhân ma trận đòi hỏi 2n3 phép 3 26
16n3 toán. Như vậy, để thực hiện phương pháp (2.6), ta cần phép toán, còn để thực 3 7n3 hiện phương pháp (1.11) ta chỉ cần phép toán. 3 Bởi vì chúng ta đã thực hiện tuyến tính hóa trong  -phương pháp, nên lược đồ (2.6) không thể có bậc hội tụ vượt quá 2. Điều này có thể được cải tiến đáng kể khi hệ (2.1) là hệ phương trình tuyến tính. 2.1.4. Thử nghiệm số Ví dụ 1.1 Hệ phương trình 1 1 2 x1 t  2 x 1 t  x 2 t , 2 x2 t x 1 t x 2 t x 2 t , t 0,1 , x1 0 1, x 2 0 1, 2t t có nghiệm chính xác là x1 e,. x 2 e Dùng phương pháp của công thức (2.6) để giải số phương trình trên. 1 Gọi e r -là chuẩn Euclid của sai số tại t 1 với bước lưới h . 80 80 Mã nguồn và Chương trình tính toán trên MATLAB được cho trong phần Phụ lục. Kết quả tính toán được cho trong bảng dưới đây.  e r80 10 5 1.7.10 5 10 6 1.7.10 5 10 7 1.7.10 5 10 8 1.7.10 5 27
2.2. Phương pháp không cổ điển giải số hệ phương trình vi phân tuyến tính cấp một Xét trường hợp khi hàm f(,)()() x t B t x g t , tức là khi (1.1) trở thành hệ phương trình vi phân tuyến tính xtBtxtgt () ()() (), t  0,1, x (0) x0 , (2.12) trong đó B t là ma trận cỡ n n, g t là hàm vectơ n chiều có các phần tử khả vi liên tục đến bậc hai (thuộc lớp C 2 0,1). 2.2.1. Phương pháp một bước 2.2.1.1. Phương pháp một bước Để giải (2.12), ta cũng bắt đầu đi từ phương pháp- và phương án một tựa của nó (hai phương pháp này đều có cấp chính xác bằng một): xi 1 xhcBx i  1 i 1 i 1 (1 cBx 1 ) i i hcg 1 i 1 (1 cg 1 ) i  và xi 1 xhBt i ( i chcx 2 ) 2 i 1 (1 cx 2 ) i hgt ( i ch 2 ) , trong đó 0 c1,c2 1. Hai công thức này có thể viết gộp lại thành một hệ phương trình đại số tuyến tính có ẩn là xi 1 : E hc B E h(1 c )B c g (1 c )g 1 i 1 x 1 i x h 1 i 1 1 i . (2.13) i 1 i E hc2 B(ti c2h) E h(1 c2 )B(ti c2h) g(ti c2h) Nhận thấy rằng hệ (2.13) nói chung không có nghiệm theo nghĩa cổ điển. Để giải (2.13) ta nhân cả hai vế của (2.13) với ma trận E hcB1i 1 c 3 E hcBt 2 i ch 2 cấp n 2n . Khi đó (2.13) trở thành: 2 2 (E hc1Bi 1) c3(E hc2B(ti c2h)) }xi 1 (E hc1Bi 1)(E h(1 c1)Bi ) c3 (E hc2B(ti c2h)(E h(1 c2 )B(ti c2h))}xi (2.14) h (E hc1Bi 1)(c1gi 1 (1 c1)gi } hc3(E hc2B(ti c2h))g(ti c2h). 28
Bổ đề 2.1 Công thức (2.14) ổn định với mọi bộ c1, c2 , c3 và có bậc hội tụ tối thiểu là một. Chứng minh Do ma trận đứng trước xi 1 là không suy biến nên (2.14) có thể viết dưới dạng xi 1 Mi xi i , trong đó 2 2 1 Mi ( EhcB 1 i 1 ) cEhcBt 3 ( 2 ( i ch 2 )) } (EhcB 1i 1 )( Eh (1)) cB 1 i cEhcBt 3 ( 2 ( i chEh 2 )( (1)( cBt 2 i ch 2 ))}, 1  (E hcB )2 cE ( hcBt ( ch )) 2 } i 1 i 1 3 2 i 2 hEhcB ( 1i 1 )( cg 1 i 1 (1)} cg 1 i hcEhcBt 3 ( 2 ( i chgt 2 ))( i ch 2 )}. Vì các phần tử của B() t và g() t thuộc lớp C 2 0,1 nên có thể chứng minh được rằng Mi 1 Lh với mọi i 0,1,2, , N . Từ đây ta có tính chất ổn định của phương pháp (2.14). Để chứng minh tính hội tụ ta nhận xét rằng, sai số địa phương của phương 2 pháp (2.14) có bậc h với mọi giá trị của c1,, c 2 c 3 . Bổ đề được chứng minh. Nhận xét Hệ (2.14) chứa ba hệ số tự do c1,, c 2 c 3 , ta tìm c1 , c2 , c3 để công thức (2.14) đạt được độ chính xác cao hơn. Viết lại hệ (2.14) như sau: 2 2 2 2 2 (1)2( cE3 hcB 1123i ccBt ( i ch 2 ))( hcB 1123 i ccBt ( i ch 2 )) x i 1 2 {(1 cEh3 )  (1 cBcB 1 )i 1 i 1 hc 1 (1 cBB 1 ) i i 1 2 2 ch3(12)( cBt 2i ch 2 ) cc 2 3 (1 chBt 2 ) ( i chx 2 )} i (2.15) h c1 gi 1 (1 c 1 ) g i c 3 g ( t i c 2 h ) 2 h c1 Bi 1( c 1 g i 1 (1 cg 1 )i ) ccBt 2 3 ( i chgt 2 ) ( i ch 2 ) . Khai triển Taylor tại điểm ti ta được: 29
h 2 c 2h 2 B B t B hB B O h 3 ; B t c h B c hB 2 B O h3 ; i 1 i 1 i i 2 i i 2 i 2 i 2 i h 2 c 2h 2 g g t g hg g O h 3 ; g t c h g c hg 2 g O h3 ; i 1 i 1 i i 2 i i 2 i 2 i 2 i h 2 h3 x x t x hx x x O h 4 . i 1 i 1 i i 2 i 6 i Thay vào (2.15) ta được: 2 2 2 h c2 h 1 cE3 2 hcBhB 1 i i B i ccBchB 2 3 i 2 i B i 2 2 h2 h 3 hcBhB2 2()()() 2 ccBchB 2 2 xhx x x 1i i 2 3 i 2 i  i i2 i 6 i 2 2 h 2 h (1 cEh3 ) (1 cBcBhB 1 )i 1 i i B i hc 1 (1 cB 1 ) i Bi hB i B i 2 2 2 2 c2 h 2 2  ch3(12)( cBchB 2i 2 i B i ) cch 2 3 (1)( cBchB 2 i 2 i )  x i 2  h2 c2 h 2 hcghg ( g ) (1 cgcgchg ) ( 2 g )] 1i i2 i 1 i 3 i 2 i 2 i 2 4 h cBhBcg1(i i )(( 1 i hg i )(1)) cg 1 i ccBchBg 2 3 ( i 2 i )( i c2 hg i )] O ( h ). hay (1)2( cE hcccBh ) 2 (22)( c ccB 2 c 2 ccB 2 )] 2 3 1 2 3i 1 2 3 i 1 2 3 i h2 h 3 h3 ( cccB 3 ) (2 c 2 2 ccBBxhx 3 ) ]}( x x } 1 3 2i 1 3 2 i i i i2 i 6 i 2 2 (1)cEh3 (12 cc 1 3 12) cBh 2 i (( cBc 1 i 1 (1) cBcc 1 i 3 2 (12) cB 2 i 2 2 3 c1 c3 c 2 c3 c 2 (1 c2))Bi h Bc i 1 (1) cBB 1 i i (12) cBcc 2 i 3 2 (1)2 ccBBx 2 2 i i } i 2 2 2 h(1 cg3 )i h ( c 1 ccg 3 2 ) i cBg 1 i i ccBg 3 2 i i  2 3 c1 c3 c 2 2 2 2 4 h (). gi cBg1 1 i cBg 1 i i ccBg 3 2 i i ccBg 3 2 i i Oh 2 2 Hệ thức trên đúng với mọi h . So sánh hệ số của h trong hai vế ta được 0 Hệ số của h : (1 c3 ) (1 c 3 ) xi 0 . 30
Hệ số của h1 : (1)2 cx3i cccBx 1 3 2 i i (12 cc 1 3 2) ccBx 3 2 i i (1)0 cg 3 i hay (1 c3 ) xi B i x i g i 0 . Hệ số của h2 : 1 c 3 x 2( cccBx ) 2( c ccBx2 ) ( c 2 ccBx 2 ) 2 2 i1 3 2 i i 1 32 i i 1 3 2 i i 2 (ccc132 1 2 cBx 2 )i i ( c 11322 1 c cc 1 cBx ) i i ( cccgcBgccBg 132 ) i 1 i i 32 i i Từ hệ phương trình x ()()()() t B t x t g t (hệ (2.1)) ta suy ra x t B t x t B t x t g t B t x t B t B t x t g t g t (2.16) B t x t B 2 t x t B t g t g t . Hệ số của h2 được viết lại như sau: 1 c 3 (Bx Bx2 Bg g ) 2( c ccBBx ) ( g ) 2 iiiiiii1 3 2 iiii ()()()cc cBx c ccBx2 c ccg cBg ccBg 321i i 132 i i 132 i 1 i i 32 i i 1 c (3 c c c ) B x B2 x B g g 0 2 1 3 2 i i i i i i i Suy ra hệ số của h2 bằng 0 khi (2.7) được thỏa mãn. Hệ số của h3 : 1 c 1 3 x 2( cccBx ) 2( cccB2 )( cccBx 2 2 2 ( cccB 3 )2( cccBBx 2 3 ) 6i132 i 2 i 132 i 132 i i 132 i 132 i i i c c c2 1 Bc (1) cBB 3 2 (12) cBcc 22 (1 cBBx ) i1 1 i i 2 i2 3 2 i i i 2 2 c c c2 ()(1 3 2 g c2 c c2)()Bg c c c 2 Bg 2 2 i 1 3 2 i i1 3 2 i i Từ (2.16) suy ra 31
x t B t x t B t x t 2B t B t x t B 2 t x t B t g t B t g t g t B t x t B t B t x t g t 2B t B t x t B 2 t B t x t g t (2.16’) B t g t B t g t g t B3 t x t B t x t 3B t B t x t 2B t g t B 2 t g t B t g t g t nên hệ số của h3 được viết lại như sau 1 c 3 (Bx3 Bx 3 BBx 2 Bg Bg 2 Bg g ) 6 ii ii iii ii ii ii i 3 2 2 (c1 cc 3 2 )( BBxiiiiiiiii Bx Bg Bg ) 2( c 1 ccBBx 3 2 ) iii 2 2 2 3 2 2 2 2(c1 ccBg 3 2 )i i ( c 1 ccBx 3 2 ) i i ( c 1 ccBg 3 2 ) i i c c ccBx3 2 c 2 2 ccBBx 3 1 B x c 1 c B B x 1 3 2 i i 1 3 2 i i i 2 i i1 1 i i c c2c c c 2 3 2 1 2c B x 2 c c2 1 c B B x (1 3 2 ) g 22i i 3 2 2 i i i 2 2 i 2 2 2 ()c1 c 3 c 2 Bi g i c 1 c 3 c 2 B i g i 1 c3 2 2 3 2 c1 c3c2 c1 c3c2 Bi xi Bi g i Bi g i 6 1 c c c c 2 1 c 3 1 3 2 B x 2B g g 3 c 2 2c c c B B x . i i i i i 1 1 3 2 i i i 6 2 2 2 Với điều kiện (2.16) thì hệ số của h3 được viết lại như sau 2c1 2 c 3 c 2 2 2 3 2 c1 cc 3 2 c 1 cc 3 2 Bxi i Bg i i Bg i i 6 2 2c1 2 c 3 c 2c1 c 3 c 2 2 c 1 2 c 3 c 2 2 BxBggi i 2 i i i c1 2 cccBBx 1 3 2 i i i 6 2 2 2 2 2 2 2 3 2 c1 cc 3 2 c 1 cc 3 2 BxBgi i i i Bg i i 3 3 1 1 1 2 2 c1 c 3 c 2 c 3c2 Bxi i 2 Bg i i g i c 1 cBBx 1 i i i 6 3 2 Hệ số của h3 bằng 0 khi 32
2 2 2 2 2c 2 c c c c2 c c2 0 c c c c 1 c 3 1 3 3 2 1 3 2 3 2 3 3 2 3 1 1 1 1 2 2 2 c1 c1 c3c2 c3c2 0 c3c2 c3c2 6 3 2 3 3 2 2 c1 c1 0 c1 c1 0 2 Suy ra nếu chọn c 0, c , c 3 thì sai số địa phương của (2.14) là bậc bốn. 1 2 3 3 Vậy (2.14) chính xác bậc ba và có dạng 2 2 2  E 3( E hB ( ti h ))  x i 1 3 3  2 2h 2  (EhB i ) 3( E hBt ( i hE )( Bt ( i hx ))  i (2.17) 3 3 3 3  2 2 2 h{ g 3( E hB ( t h )) g ( t h )}. i3 i 3 i 3 2.2.1.2. Phương trình thử Áp dụng công thức (2.17) cho phương trình thử: x (t)  x(t), t 0,1, x(0) x0 . Ta có B(); t   B t x t g t  x t g( t ) 0. Khi đó (2.17) được viết lại thành: 22 2 1 EEhx 3(  ) i 1 ( Eh  ) 3( EhEhx  )(  ) i . 3 3 3 Suy ra 4 2 2 2 4 (h) 4(h) xi 1 4 (h) xi 3 3 hay với z h ta có xi 1 R z xi , trong đó 33
1 z 2 / 6 R(z) . 1 z z 2 /3 Kiểm tra điều kiện R(z) 1 ta thấy R z chứa cả nửa trái mặt phẳng phức. Điều này có nghĩa là công thức (2.14) là phương pháp một bước, một nấc, có chứa phương pháp là ổn định -A và chính xác bậc ba. Nhận xét rằng, nếu các tham số c1,, c 2 c 3 được chọn thỏa mãn điều kiện (2.7) thì ta có lược đồ sai phân chính xác cấp hai. 1 Đặc biệt, chọn c ,c 0 thì (2.14) trở thành: 1 2 3 1 1 1 1 1 1 (E hB )2 x (E hB )(E hB )x h(E B )( g g ) 2 i 1 i 1 2 i 1 2 i i 2 i 1 2 i 1 2 i 1 x x h(B x B x g g ) . i 1 i 2 i 1 i 1 i i i i Đây chính là phương pháp hình thang. Có thể kiểm tra được rằng, phương pháp (2.14) cho phương trình thử và phương pháp Runge-Kutta với bảng Butcher 0 ½ -1/2 2/3 1/6 ½ ¼ ¾ cho cùng một kết quả. Tuy nhiên, với B() t là ma trận phụ thuộc vào t thì kết quả tính toán theo phương pháp Runge-Kutta với bảng trên và theo phương pháp (2.14) là khác nhau. Nghĩa là chúng là những phương pháp khác biệt. Nhận xét 7n3 Để thực hiện phương pháp (2.3) cần phép toán số học ( 2n3 phép toán cần cho 3 n3 nhân ma trận và phép toán giải hệ phương trình tuyến tính). Để thực hiện 3 8n3 2n3 phương pháp Runge-Kutta hai nấc cần phép toán và phép toán cho 3 3 34
phương pháp đường chéo hiển. Như vậy, trong trường hợp tổng quát, phương pháp hai nấc   0 1  1 2  1/2 ½ 1 3 với  có cấp chính xác bậc ba, đòi hỏi số phép toán số học ít hơn phương 2 6 pháp (2.14). Tuy nhiên, nếu có thể tính ma trận B2 () t một cách giải tích (dễ dàng làm được khi các phần tử của B() t là các đa thức và ghi vào ô nhớ, thì phương pháp (2.14) có lợi thế hơn so với phương pháp Runge-Kutta hai nấc: Phương pháp (2.14) n3 chỉ cần phép toán số học. 3 Thử nghiệm số Ví dụ 2.1 Phương trình 501 500 0 x t x t , x 0 , 500 501 2 t 1001t t 1001t có nghiệm chính xác là x1 e e x2 e e . Dùng phương pháp của (2.6) để giải số phương trình trên. Gọi e rr -là chuẩn Euclid của sai số tại t 1. Các tính toán được thực hiện trên MATLAB. Mã nguồn và Chương trình được cho trong phần Phụ lục. Kết quả tính toán được cho trong bảng dưới đây. h e rr 0.1 0.1.10 2 0.05 0.26.10 5 0.025 0.33.10 6 2.2.2. Phương pháp đa bước 35
Mục này trình bày cách xây dựng các lược đồ sai phân đa bước không cổ điển giải số hệ (2.12) do Bulatov đề xuất (xem [11]). Phương pháp này cho phép xây dựng các lược đồ sai phân với ưu thể hơn các phương pháp cổ điển. Thí dụ, từ phương pháp này có thể xây dựng được các lược đồ sai phân cấp ba hoặc cấp bốn, tương ứng hai hoặc ba bước là ổn định - A. Để giải bài toán (2.12) theo phương pháp đa bước ta bắt đầu đi từ m phương pháp khác nhau có cùng k bước theo công thức sau: k k l l  jx i 1 j h   j B( t i 1 j ) x i 1 j g ( t i 1 j ), l 1,2, , m . (2.18) j 0 j 0 Ở đây chỉ số l là chỉ số của phương pháp. Trong các lược đồ (2.18) có thể có các lược đồ ổn định cũng có thể có các lược đồ không ổn định. Viết lại (2.18) thành hệ phương trình đại số tuyến tính có ẩn là xi 1 như sau: 1 1 1 1 1 E h B j E h j Bi 1 j  j g i 1 j 0 0 i 1 2 2 k 2 2 k 2 0 E h 0 Bi 1 j E h j Bi 2 j  j g i 1 j xi 1  xi 1 j h . (2.19) j 1 j 0 m E h m B m m m 0 0 i 1 j E h j Bi 1 j  j g i 1 j Hệ (2.19) nói chung không có nghiệm theo nghĩa cổ điển. Để giải (2.19) ta nhân cả hai vế của nó với ma trận: 1 1 Li 1 ( 0 h  0 B i 1 c 1 E c m 1 E ) (2.20) cấp n mn , trong đó c1, c 2 , , cm 1 là những hằng số tùy ý. Sau khi nhân với ma trận Li 1 , hệ (2.19) trở thành: k Mi 1 x i 1  M i 1 j x i 1 j  i 1 , (2.21) j 1 trong đó 36
k 1  j gi 1 j 1 1 1 1 j 0 E h B j E h j Bi 1 j 0 0 i 1 k 2 2 2 2  2 g 0 E h0 Bi 1 j E h j Bi 2 j  j i 1 j Mi 1 Li 1 , M i 1 j Li 1 , i 1 hLi 1 j 0 . mE h mB m m 0 0 i 1 j E h j Bi 1 j k  m g  j i 1 j j 0 Vì mỗi lược đồ trong (2.19) có cấp chính xác l nên lược đồ (2.21)cũng có cấp chính xác l với mọi c1, c 2 , , cm 1 . Phương pháp chung để xây dựng và khảo sát công thức (2.21) dùng để giải số bài toán (2.12) là: 1. Bắt đầu đi từ công thức (2.18). 2. Xây dựng công thức dạng (2.21) 3. Phân tích (2.21) thành chuỗi Taylor và chọn các hệ số c1, c 2 , , cm 1 sao cho lược đồ (2.21) có bậc l p , trong đó p m 1, ngoài ra còn r m 1 p tham số tự do. 4. Chọn r tham số còn lại để được lược đồ ổn định và trong trường hợp lý tưởng, là ổn định - A. Nhận xét Nếu nhân hai vế của hệ (2.19) với ma trận Li 1 (d1E d 2E d m E) (2.20a) cấp n mn, trong đó d1,d2 , ,dm là các hằng số tùy ý, thì một lần nữa ta lại nhận được một lược đồ sai phân k bước cổ điển có dạng: k k  j xi 1 j h j f (xi 1 j ,ti 1 j ). j 0 j 0 Khác biệt cơ bản của phương pháp do Bulatov đề xuất là Li 1 xác định bởi công thức (2.20) phụ thuộc vào h và Bi 1 , trong khi đó phương pháp cổ điển (2.20a) không phụ thuộc vào h và Bi 1 . Dưới đây sẽ minh họa phương pháp nêu trên qua các lược đồ cụ thể. 37
Lược đồ 2.1 Xây dựng lược đồ hai bước chính xác bậc ba giải bài toán (2.12) Để giải phương trình (2.12) bằng phương pháp hai bước, trước hết ta xây dựng một lược đồ sai phân hai bước bậc hai như sau. Theo khai triển Taylor tại ti 1 ta có: 4h2 x x t x t 2h x 2hx x O(h3 ) ; i 1 i 1 i 1 i 1 i 1 2 i 1 h 2 x x t x t h x hx x O(h 3 ). i i i 1 i 1 i 1 2 i 1 Suy ra 3 xi 1 4xi 3xi 1 2hxi 1 O(h ) hay 3 3xi 1 4 x i x i 1 2 h ( B i 1 x i 1 g i 1 ) O ( h ) a Theo khai triển Taylor tại ti ta có: 2 h 3 xi 1 x ti 1 x t i h xi hxi xi O(h ), 2 h 2 x x t x t h x hx x O(h 3 ). i 1 i 1 i i i 2 i Suy ra x x 2hx O(h 3 ) i 1 i 1 i 3 xi 1 xi 1 2h(Bi xi g i ) O(h ). b Theo khai triển Taylor tại ti 1 ta có: h2 x x t x t h x hx x O(h3); i i i 1 i 1 i 1 2 i 1 4h2 x x t x t 2h x 2hx x O(h3 ). i 1 i 1 i 1 i 1 i 1 2 i 1 Suy ra 4x x 3x 2hx O(h3 ) i i 1 i 1 i 1 3 xi 1 4xi 3xi 1 2h(Bi 1 xi 1 gi 1 ) O(h ). c 38
Từ (a), (b), (c) ta có một lược đồ sai phân hai bước cấp hai như sau: 3xi 1 4xi xi 1 2h(Bi 1xi 1 gi 1); xi 1 xi 1 2h(Bi xi gi ); (2.22) xi 1 4xi 3xi 1 2h(Bi 1xi 1 gi 1). Đưa công thức (2.22) về dạng công thức (2.21) (3E 2hBi 1)xi 1 4xi xi 1 2hgi 1 xi 1 2hBi xi xi 1 2hgi xi 1 4xi (3E 2hBi 1)xi 1 2hgi 1 hay 3E 2hB 4E E 2hg i 1 i 1 E xi 1 2hBi xi E xi 1 2hgi . (2.22a) E 4E 3E 2hBi 1 2hgi 1 Nhân hai vế của hệ phương trình trên với ma trận (3E 2hBi 1 c1E c2E) ta được: 2 (3E 2hBi 1 ) (c1 c2 )E)]xi 1 (12 4c2 )E 2h(c1Bi 4Bi 1 )xi (2.23) ( 3 c1 3c2 )E 2h(Bi 1 c2 Bi 1 )xi 1 2h(3E 2hBi 1 )g i 1 c1g i c2 g i 1 ] hay 2 2  9 c1 c2 E 12hBi 1 4h Bi 1 xi 1  12 4c2 E 2h c1Bi 4Bi 1 xi (2.23’)  3 c1 3c2 E 2h Bi 1 c2 Bi 1 xi 1 2h 3E 2hBi 1 gi 1 c1gi c2 gi 1 . Khai triển Taylor tại ti 1 ta được: h2 h3 x x t x hx x x O(h4), i i i 1 i 1 2 i 1 6 i 1 4h2 8h3 x x t x 2hx x x O(h4), i 1 i 1 i 1 i 1 2 i 1 6 i 1 h2 4h2 B B t B hB B O(h3);B B t B 2hB B O(h3), i i i 1 i 1 2 i 1 i 1 i 1 i 1 i 1 2 i 1 h2 4h2 g g t g hB g O(h3);g g t g 2hg g O(h3). i i i 1 i 1 2 i 1 i 1 i 1 i 1 i 1 2 i 1 Thay vào (2.23’) ta được 39
2 4h 2 4h 2 (9 c c )E 12h B 2hB B 4h 2 B 2hB B 1 2 i 1 i 1 i 1 i 1 i 1 i 1 2 2 4h 2 8h3 xi 1 2hxi 1 xi 1 xi 1 2 6 h2 4h2  (12 4c2 )E 2h c1 (Bi 1 hBi 1 Bi ) 4(Bi 1 2hBi 1 Bi 1 )  2 2  h2 h3 xi 1 hxi 1 xi xi 1 2 6 4h2  ( 3 c1 3c2 )E 2h Bi 1 2hBi 1 Bi 1 c2 Bi 1 xi 1 2  4h 2 4h2 4h 2  3(g i 1 2hgi 1 g i 1 ) 2h Bi 1 2hBi 1 Bi 1 gi 1 2hgi 1 g i 1 2 2 2 4 . 2h  O h h 2 c g hg g c g 1 i 1 i 1 i 1 2 i 1 2  Chú ý rằng ta chỉ cần đến khai triển của h3 nên các biểu thức chứa h4 được đưa vào công thức O h 4 do đó không có trong công thức dưới đây, khi ấy công thức (2.23’) được viết lại như sau 2 2 3  9 c1 c2 E h 12Bi 1 h 24Bi 1 4Bi 1 h 24Bi 1 16Bi 1Bi 1  4h2 8h3 xi 1 2hxi 1 xi 1 xi 1 2 6 h2 h3 2 3 (12 4c2 )E h 2c1 8 Bi 1 h (2c1 16)Bi 1 h (c1 16)Bi 1 xi 1 hxi 1 xi 1 xi 1 2 6 2 3  3 c1 3c2 E h 2 2c2 Bi 1 h .4Bi 1 h 4Bi 1xi 1 2h 3 c1 c2 gi 1 2 3 4 h 12gi 1 4Bi 1gi 1 2c1gi 1 h 12gi 1 8Bi 1gi 1 8Bi 1gi 1 c1gi 1 O h 0 Hệ số của h : (9 c1 c2 )xi 1 (12 4c2 )xi 1 ( 3 c1 3c2 )xi 1 0 . Hệ thức này luôn đúng với mọi c1, c 2 . Hệ số của h1 : (9 c1 c2 ).2xi 1 12Bi 1 xi 1 (12 4c2 )xi 1 (2c1 8)Bi 1xi 1 (2 2c2 )Bi 1xi 1 6 2c1 2c2 (6 2c 2c )x (6 2c 2c )B x (16 2c 2c )g 1 2 i 1 1 2 i 1 i 1 1 2 i 1 (6 2c1 2c2 )xi 1 (6 2c1 2c2 )xi 1 0 40
Hệ thức này luôn đúng với mọi c1, c 2 . Hệ số của h2 : 1 (9 c c ).2x 12B .2x ( 24B 4B2 )x 12 4c . x 2c 8 B x 1 2 i 1 i 1 i 1 i 1 i 1 i 1 2 2 i 1 1 i 1 i 1 2c1 16 Bi 1xi 1 4Bi 1xi 1 12 2c1 gi 1 4Bi 1gi 1 2 12 2c1 xi 1 12 2c1 Bi 1xi 1 16 2c1 Bi 1 xi 1 12 2c1 gi 1 4Bi 1gi 1 4Bi 1 xi 1 2 12 2c1 xi 1 12 2c1 Bi 1xi 1 16 2c1 Bi 1 Bi 1xi 1 gi 1 12 2c1 gi 1 4Bi 1gi 1 4Bi 1xi 1 2 12 2c1 xi 1 12 2c1 Bi 1xi 1 Bi 1xi 1 Bi 1gi 1 gi 1 Do (2.16) nên hệ số của h2 là 12 2c1 xi 1 12 2c1 xi 1 0 . Hệ thức này luôn đúng với mọi c1, c 2 . Hệ số của h3 : 8 (9 ccx ). ( 12 Bx ).2 ( 24 BBx 42 ).2 24 Bx 16 BBx 1216 i ii 11 iii 111 ii 11 iii 111 1 1 124 cx 28 cBx 216 c Bxc 16 BxBx 4 2116i i 111 2 i i 111 i i 11 i i 11 i 12 c1 gi 1 8 BB i 1 g i 1 8 Bi 1g i 1 4 2 2 10ccx1211 i c 20 Bx i 11 i 2 c 1 32 BxBx i 11111 i 8 i i c 12 Bx i 11 i 3 3 12 c1 gi 1 8 BB i 1 g i 1 8 B i 1 g i 1 16 B i 1 B i 1 x i 1 Do (2.16) nên hệ số của h3 bằng 4 2 2 10 ccx1211 i c 20 BBxBxBgg iiiiiiii 11111111 2 c 1 32 BBxg iiii 1111 3 3 2 8BBxgiiii 1111 c 1 12 Bx ii 11 12 cg 11 iiiii 8 Bg 11 8 Bg 11 16 BBx iii 111 4 2 3 2 10 c1 c 2 xi 1 12 cBxBx 11111 iiii 3 BBx iii 111 2 BgBgBgg iiiiiii 1111111 3 3 Do (2.16’) nên suy ra hệ số của h3 bằng 0 khi 4 2 10 c c 12 c c 2 c 6 . (2.24) 31 3 2 1 1 2 Với điều kiện (2.24) thì lược đồ (2.23) là họ lược đồ sai phân hai bước chính xác bậc ba, tuy nhiên nó chứa cả các phương pháp ổn định và các phương pháp không ổn định. 41
Áp dụng lược đồ (2.23) cho phương trình ổn định x 0 ta được: 9 c c x 12 4c x 3 c 3c x 1 2 i 1 2 i 1 2 i 1 9 c1 c2 xi 1 12 4c2 xi 3 c1 3c2 xi 1 0 Phương trình sai phân này có phương trình đặc trưng tương ứng là 2 9 c1 c2  12 4c2  3 c1 3c2 0 3 c1 3c2 Phương trình này có nghiệm 1 1; 2 . 9 c1 c2 Lược đồ (2.23) ổn định khi phương trình đặc trưng có nghiệm thoả mãn  1 tức là 3 c 3c 1 2 1. Bất phương trình này tương đương với 9 c1 c2 c2  3 c1 ,3. (2.25) Đặc biệt, nếu chọn c1 6, c2 0 thì lược đồ (2.23) được viết lại như sau: 3E 2hB 2 6E x  i 1  i 1 ( 2.26) 12E 2h(3Bi 2Bi 1 )xi 3E 2hBi 1 xi 1 2h(3E 2hBi 1 )gi 1 6g i ]. Áp dụng lược đồ (2.26) cho phương trình thử x (t)  x(t), t 0,1, x(0) x0 ta được 2  3E 2h 6Exi 1 12E 2hxi 3E 2hxi 1 0 . Phương trình đặc trưng tương ứng là 15 6z z 2  2 12 2z  3 z 0 với z 2h . Hình 2.1 Miền ổn định của hệ sai phân với z 2h được mô tả trên Hình 2.1 (miền gạch là miền không ổn định. Như vậy phương pháp này là ổn định - A. Nhận xét 1 4 1 Nếu nhân hai vế của (2.22a) với ma trận L ( E E E) thì ta thu được phương i 1 6 6 6 pháp Milne 42
h x x (B x 4B x B x g 4g g ) . i 1 i 1 3 i 1 i 1 i i i 1 i 1 i 1 i i 1 Lược đồ 2. 2 Trong lược đồ này ta xây dựng một phương pháp ba bước chính xác bậc bốn giải bài toán (2.12). Đầu tiên ta xây dựng một lược đồ sai phân ba bước chính xác bậc ba như sau. Theo khai triển Taylor tại ti 1 ta có: h 2 h3 x x t x t h x hx x x O(h 4 ) ; i i i 1 i 1 i 1 2 i 1 6 i 1 4h 2 8h3 x x t x t 2h x 2hx x x O(h4 ) ; i 1 i 1 i 1 i 1 i 1 2 i 1 6 i 1 9h2 27h3 x x t x t 3h x 3hx x x O(h4 ) . i 2 i 2 i 1 i 1 i 1 2 i 1 6 i 1 Suy ra 4 18xi 9xi 1 2xi 2 11xi 1 6hxi 1 0(h ). (a) Khai triển Taylor tại ti ta có: h 2 h 3 x x t x t h x hx x x O(h 4 ) ; i 1 i 1 i i i 2 i 6 i h 2 h3 x x t x t h x hx x x O(h 4 ) ; i 1 i 1 i i i 2 i 6 i 4h 2 8h3 x x t x t 2 x 2hx x x O(h4 ) . i 2 i 2 i i i 2 i 6 i Suy ra 4 2xi 1 6xi 1 xi 2 3xi 6hxi O(h ) . (b) Theo khai triển Taylor tại ti 1 ta có: 4h 2 8h 3 x x t x t 2h x 2hx x x O(h 4 ) ; i 1 i 1 i 1 i 1 i 1 2 i 1 6 i 1 h2 h3 x x t x t h x hx x x O(h4 ) ; i i i 1 i 1 i 1 2 i 1 6 i 1 43
h 2 h 3 x x t x t h x hx x x O(h 4 ) . i 2 i 2 i 1 i 1 i 1 2 i 1 6 i 1 Suy ra 4 xi 1 6xi 2xi 2 3xi 1 6hxi 1 O(h ) . (c) Theo khai triển Taylor tại ti 2 ta có: 9h 2 27h3 x x t x t 3h x 3hx x x O(h 4 ) ; i 1 i 1 i 2 i 2 i 2 2 i 2 6 i 2 4h 2 8h3 x x t x t 2h x 2hx x x O(h 4 ) ; i i i 2 i 2 i 2 2 i 2 6 i 2 h 2 h 3 x x t x t h x hx x x O(h 4 ) . i 1 i 1 i 2 i 2 i 2 2 i 2 6 i 2 Suy ra 4 2xi 1 9xi 18xi 1 11xi 2 6hxi 2 O(h ) . (d) Từ (a), ( b), (c), (d) ta có một lược đồ sai phân ba bước chính xác bậc ba như sau: 11xi 1 18 x i 9 x i 1 2 x i 2 6 h B i 1 x i 1 g i 1 2xi 1 3 x i 6 x i 1 x i 2 6 h B i x i g i (2.27) xi 1 6 x i 3 x i 1 2 x i 2 6 h B i 1 x i 1 g i 1 2xi 1 9 x i 18 x i 1 11 x i 2 6 h B i 2 x i 2 g i 2 Đưa (2.27) về dạng (2.21) (11E 6 hBi 1 ) x i 1 18 x i 9 x i 1 2 x i 2 6 hg i 1 2xi 1 ( 3 E 6 hB i ) x i 6 x i 1 x i 2 6 hg i xi 1 6 x i (3 E 6 hB i 1 ) x i 1 2 x i 2 6 hg i 1 2xi 1 9 x i 18 x i 1 (11 E 6 hB i 2 ) x i 2 6 hg i 2 hay 11E 6hB 18 9 2 g i 1 i 1 2 3E 6hBi 6 1 gi x x x x 6h 1 i 1 6 i 3E 6hB i 1 2 i 2 g i 1 i 1 2 9 18 11E 6hBi 2 gi 2 Nhân hai vế của hệ phương trình trên với ma trận: 44
Li 1 (11E 6hBi 1 c1E c2 E c3 E ) ta được 2 (11E 6hBi 1 ) (2c1 c2 2c3 )Exi 1 (198 3c1 6c2 9c3 )E 6h(c1Bi 18Bi 1 )xi ( 99 6c1 3c2 18c3 )E 6h(9Bi 1 c2 Bi 1)xi 1 (22 c1 2c2 11c3 )E 6h( c3 Bi 2 2Bi 1 )xi 2 (2.28) 6h(11E 6hBi 1 )gi 1 c1gi c2 gi 1 c3 gi 2  Theo khai triển Taylor tại ti 2 ta có: 9h2 27 h 3 81 h 4 xxtxthxhx 3 3 x x xOh 4 (5 ); i 1 i 1 i 2 i 2 i 22 i 2 6 i 2 24 i 2 4h2 8 h 3 16 h 4 xxtxt 2 hx 2 hx x x xOh 4 (5 ); i i i 2 i 2 i 22 i 2 6 i 2 24 i 2 h2 h 3 h 4 xxt xthxhx x x x 4 O( h5 ); i 1 i 1 i 2 i 2 i 22 i 2 6 i 2 24 i 2 9h2 27 h 3 BBt Bt 3 hB 3 hB B BOh (4 ); i 1 i 1 i 2 i 2 i 22 i 2 6 i 2 4h2 8 h 3 BBtBt 2 hB 2 hB B BOh (4 ); i i i 2 i 2 i 22 i 2 6 i 2 h2 h 3 BBt BthBhB B BOh (4 ); i 1 i 1 i 2 i 2 i 22 i 2 6 i 2 9h2 27 h 3 g g t g t 3 h g 3 hg g g O ( h4 ); i 1 i 1i 2 i 2 i 22 i 2 6 i 2 4h2 8 h 3 ggtgt 2 hg 2 hg g gOh (4 ); i i i 2 i 2 i 22 i 2 6 i 2 h2 h 3 ggt gthghg g gOh (4 ). i 1 i 1 i 2 i 2 i 22 i 2 6 i 2 Thay vào (2.28) ta được 45
2 2 3 ((1212 cccEhBh1 2 2) 3 132i 2 396 B i 2 36 Bh i 2 594 B i 2 216 BB i 2 i 2 4 2 h 594 Bi 2 324 B i 2 324 B i 2 B i 2 ) 2 3 9h 27 h 81 (4) xi 2 3 hx i 2 x i 2 x i 2 x i 2 2 6 24 (198 3cccEhc 6 9 ) (6 108) Bhc 2 (12 324) Bhc 3 (12 486) B 1 2 3 1i 2 1 i 2 1 i 1 h4 (8 c 81) B 1i 2 2 3 4 4h 8 h 16 h (4) xi 2 2 hx i 2 x i 2 x i 2 x i 2 2 6 24 2 (996 c1 3 c 2 18) c 3 E h (546) c 2 Bi 2 h (1626) c 2 B i 2 h3(2433) c B h 4 (243 c ) B 2i 2 2 i 2 2 3 4 h h h 4 xi 2 hx i 2 x i 2 x i 2 x i 2 2 6 24 2 3 4 22cccEhc1 2 2 11 3 612 3 Bhi 2 (36 Bh i 2 ) (54 Bh i 2 ) (54 Bx i 2 ) i 2 6h (11 c c c ) g h2 36 B g 198 g 1 2 3i 2 i 2 i 2 i 2 3 4 h 108 Bgi 2 i 2 108 Bg i 1 i 2 297 gh i 2 162 Bg i 2 i 2 162 Bg i 2 i 2 324 B i 2gi 2 297 g i 2 2 3 4 5 h 12 ccg1 6 2 i 2 h 12 ccg 1 3 2 i 2 hccg 8 1 2 i 2 Oh Hệ số của h0 : (121 2c c 2c )x 198 3c 6c 9c x 1 2 3 i 2 1 2 3 i 2 99 6c1 3c2 18c3 xi 2 22 c1 2c2 11c3 xi 2 0. Hệ thức này luôn đúng với mọi c1, c 2 . Hệ số của h1 : (121 2c1 c2 2c3 ).3xi 2 132Bi 2 xi 2 (198 3c1 6c2 9c3 ).2xi 2 (6c1 108)Bi 2 xi 2 ( 99 6c1 3c2 18c3 ).xi 2 (54 6c2 )Bi 2 xi 2 ( 6c3 12)Bi 2 xi 2 (66 6c1 6c2 6c3 )gi 2 (66 6c 6c 6c ).x (66 6c 6c 6c )B x (66 6c 6c 6c )g 1 2 3 i 2 1 2 3 i 2 i 2 1 2 3 i 2 (66 6c1 6c2 6c3 ).xi 2 (66 6c1 6c2 6c3 )xi 2 0. Hệ thức này luôn đúng với mọi c1, c 2 . Hệ số của h 2 : 46
9 (121 2c c 2c ) x 132B .3x 396B 36B2 x 1 2 3 2 i 2 i 2 i 2 i 2 i 2 i 2 (198 3c1 6c2 9c3 )2xi 2 (6c1 108).2xi 2 Bi 2 (12c1 324)Bi 2 xi 2 1 ( 99 6c 3c 18c ). x '' (54 6c )B .x (162 6c )B x 1 2 3 2 i 2 2 i 2 i 2 2 i 2 i 2 36Bi 2 xi 2 36Bi 2 gi 2 198 12c1 6c2 g i 2 2 (198 12c1 6c2 )xi 2 (198 12c1 6c2 ) Bi 2 xi 2 Bi 2 xi 2 Bi 2 gi 2 g i 2 Do (2.16) nên hệ số của h2 bằng 0. Hệ số của h3 : 27 9 (1212 c c 2) c x (132) B . x 396 B 36 B2 .3 x 1 2 36i 2 i 2 2 i 2 i 2 i 2 i 2 8 594B 216 B B x (1983 c 6 c 9). c x i 2 i 2 i 2 i 2 1 2 36 i 2 (6c1 108) Bi 2 2. x i 2 (12 c 1 324) B i 2 .2 x i 2 (12 c 1 486) B i 2 x i 2 1 1 ( 99 6c 3 c 18 c ). x (54 6 c ) B x (162 6 c ) B x 1 2 36i 2 2 i 2 2 i 2 2 i 2 i 2 243 3cBx222 i i 54 Bx i 22 i 108 Bg i 22 i 297 g i 2 12 ccg 122 3 i 108Bi 2 g i 2 297 12c1 3 c 2 xi 2 405 12 c 1 3 c 2 B i 2 x i 2 702 24c1 6 c 2 Bi 1 x i 2 2 108Bxii 12 297 12 ccBx 1222 3 ii 216 BBx iii 222 108 BgBg iiii 2222 297 12c1 3 c 2 gi 2 Do (2.16) và (2.16’) nên hệ số của h3 được viết lại như sau 2 297 12c1 3c2 xi 2 405 12c1 3c2 Bi 2 Bi 2 xi 2 Bi 2 xi 2 Bi 2 gi 2 gi 2 2 702 24c1 6c2 Bi 1 Bi 2 xi 2 gi 2 108Bi 1 Bi 2 xi 2 gi 2 297 12c1 3c2 Bi 2 xi 2 216Bi 2 Bi 2 xi 2 108 Bi 2 gi 2 Bi 2 gi 2 297 12c1 3c2 gi 2 3B B x B3 x 2B g 297 12c 3c x 297 12c 3c i 2 i 2 i 2 i 2 i 2 i 2 i 2 1 2 i 2 1 2 2 Bi 2 gi 2 Bi 2 gi 2 297 12c1 3c2 xi 2 297 12c1 3c2 xi 2 0 Hệ số của h4 : 47
81 4 272 9 (1212 c1 c 2 2) c 3 xi 1 132 B i 2 . x i 2 396 B i 2 36 B i 2 . x i 2 24 6 2 2 594Bi 2 216 B i 2 B i 1 .3 x i 2 594 B i 2 324 B i 2 324 B i 2 B i 2 x i 2 16 8 (1983 c 6 c 9). c x 4 (6 c 108) B x (12 c 324) B 2 x 1 2 324i 2 1 i 2 6 i 2 1 i 2 i 2 1 (12c 486) B .2 x (8 c 486) B x ( 99 6 c 3 c 18 c ). x 4 1i 2 i 2 1 i 2 i 2 1 2 324 i 1 1 1 (546) c B x (1626) c B x (243 3c ) B x (243 c ) B x 2i 26 i 2 2 i 2 2 i 2 2i 2 i 2 2 i 2 i 2 54Bx 162 BgBg 324 Bg 297 8 ccg i 2i 2 i 2 i 2 i 2 i 2 i 2 i 2 1 2 i 2 561 17 1 3 4 cccx1232 i 4598 ccBx 1222 i i 121524 ccBx 1222 3 i i 2 2 2 2 2 162Bxi 22 i 105324 ccBx 1222 3 i i 648 BBx i 222 i i 2978 ccBx 1222 i i 2 324Bi 2 x i 2 324 B i 2 B i 2 x i 2 162 Bi 2 g i 2 B i 2 g i 2 324 B i 2 g i 2 297 8 c 1 c 2 g i 2 Do (2.16) và (2.16’) nên hệ số của h4 được viết lại như sau 3 561 17 1 3 4 Bi 2 x i 2 3 B i 2 B i 2 x i 2 B i 2 x i 2 c1 c 2 c 3 xi 2 459 8 c 1 c 2 B i 2 2 2 2 2 2 2Bi 2 g i 2 B i 2 g i 2 B i 2 g i 2 g i 2 2 1215 24c1 3 c 2 Bi 2 B i 2 x i 2 B i 2 x i 2 B i 2 g i 2 g i 2 2 2 162Bi 2 B i 2 xi 222222 BxBgg i i i i i 1053 24 ccBBxg 122222 3 i i i i 2 648BBxi 222 i i 297 8 ccBx 1222 i i 324 Bx i 22 i 324 BBx i 222 i i 162 Bi 2 g i 2 B i 2 g i 2 324 B i 2 g i 2 297 8 c 1 c 2 g i 2 561 17 1 3 4 c1 c 2 c 3 xi 2 2 2 2 2 Bx 4 BBx 3 Bg 3 Bx 2 6 BBx 2 ii 22 iii 222 ii 22 ii 22 iii 222 4 3 2 297 8cc12 5 BBgi 222 i i Bx i 22 i Bg i 22 i 3 BgBg i 22 i 22 i B g g i 2 i 2 i 2 Từ (2.16’) ta có thể viết tiếp biểu thức trên như sau: 48
x 4 t 3B 2 t B t x t B3 t x t B t x t B t x t 3B t B t x t 3B 2 t x t 3B t B t x t 2B t g t 2B t g t 2B t B t g t B 2 t g t B t g t B t g t g t 3B 2 t B t x t B 3 t B t x t g t B t x t B t B t x t g t 3B t B t x t 3B 2 t x t 3B t B t B t x t g t 2B t g t 2B t g t 2B t B t g t B 2 t g t B t g t B t g t g t B t x t 4B t B t x t 3B t g t 3B 3 t x t 6B t B 2 t x t 5B t B t g t B 4 t x t B3 t g t 3B t g t B 2 t g t B t g t g t Hệ số của h4 là 561 17 1 3 4 4 c1 c2 c3 xi 2 297 8c1 c2 xi 2 . 2 2 2 2 Suy ra hệ số của h4 bằng 0 khi 561 17 1 3 c1 c 2 c 3 297 8 c 1 c 2 0 c 1 c 2 3 c 3 33 . (2.29) 2 2 2 2 Với điều kiện (2.29) thì lược đồ (2.28) là lược đồ ba bước chính xác bậcbốn, tuy nhiên nó chứa cả những lược đồ ổn định và các lược đồ không ổn định. Áp dụng lược đồ (2.29) với c1 33, c 2 c 3 0 cho phương trình thử x () t  x (), t t  0,1, x (0) x0 , ta được 116Eh  2 2 ccc 2  (996318)6(9 cccEhcx   ) 1 2 3 1 2 3 2i 1 (22 c1 2 c 2 11) c 3 E 6( h c 3 2)   xi 2 . Phương trình đặc trưng tương ứng là 18722 z z2  3 9915 z  2 999 z  112 z 0 với z 6 h . Imz Miền ổn định của phương pháp được mô tả trên Hình 2.2. (miền gạch là miền không ổn định). 49
Do đó phương pháp là ổn định - A, hai hoặc ba bước tương ứng bậc ba hoặc bậc bốn. Hình 2.2 50
CHƯƠNG 3 Về một phương pháp không cổ điển giải số hệ phương trình vi phân cấp hai Tương tự như trong trường hợp hệ hệ phương trình vi phân cấp một đã trình bày trong Chương 2, trong Chương này chúng tôi trình bày một phương pháp không cổ điển do M. V. Bulatov và G. V. Berghe đề xuất giải hệ phương trình vi phân cấp hai (xem [4], 2009). Một số phương pháp giải hệ phương trình vi phân cấp hai đã được trình bày trong [5], Chương 3.10 (trang 461-474); trong [7], Chương 2.12 (trang 114-120; v.v., Hiện nay phương pháp đa bước cũng như các biến thể một cột (one leg versions) của nó đã được nghiên cứu nhiều và đã được áp dụng để giải hệ phương trình vi phân cấp hai tuyến tính và phi tuyến. Phương pháp k bước bậc cao thường được xây dựng dưới dạng tổ hợp tuyến tín của một số phương pháp k bước bậc thấp hơn. Nội dung cơ bản của phương pháp do Bulatov đề xuất dựa trên việc xây dựng ma trận tổ hợp theo ma trận ban đầu A() t , bước lưới h và các tham số. Tổ hợp này cho phép xây dựng một họ các lược đồ sai phân hai bước ổn định - P bậc bốn. 3.1. Phương pháp không cổ điển giải số hệ phương trình vi phân tuyến tính cấp hai 3.1.1. Phương pháp cổ điển Cho phương trình x t A t x t f t , t 0,1 3.1 x 0 x0 , x 0 x0 , trong đó A t là ma trận cấp n n, f t là một hàm vectơ n chiều. Các phần tử của ma trận A t và vectơ f t là những hàm số đủ trơn. 51
Xây dựng lưới điểm đều ti ih; i 0,1, , N ; Nh 1 trên đoạn [0,1]. Kí hiệu xi là giá trị xấp xỉ của x ti ,, A i A t i f i f t i . Phương pháp tuyến tính k bước cổ điển có dạng: k k 2  jx i 1 j h   j A i 1 j x i 1 j f i 1 j . j 0 j 0 Nếu j k j và j  k j ,j 0,1, , k /2 thì công thức trên được gọi là công thức đối xứng. Chúng ta đã biết rằng công thức k bước có cấp chính xác không vượt quá k 2 nếu k là số chẵn và không vượt quá k 1 nếu k là số lẻ (xem [5]. Định lý 10.4, trang 468). Khi k 2 công thức trên chính là công thức Numerov chính xác bậc bốn: h2 x 2 xx Axf 10 AfAxf . i 1 ii 112 iii 1 1 1 iiiii 1 1 1 3.1.2. Lược đồ sai phân mới Ký hiệu các toán tử vi phân xi  xi 1 2xi xi 1 , j j j  j xi   0 xi 1 1 xi  2 xi 1 , j 1, 2 . Xét lược đồ sai phân hai bước giải bài toán (1.1): 2 xi  h 1 Ai xi f i  . hay 2 1 1 1 xiii 1 2 xxh 1  01111 Axf iii  Axf iii  2111 Axf iii (3.2) Khai triển Taylor tại điểm ti ta được: 2 h 3 xi 1 x t i 1 x t i h xi hxi xi O h ; 2 h 2 x x t x t h x hx x O h 3 . i 1 i 1 i i i 2 i Theo (1.1) ta có 52
x x t A x f i 1 i 1 i 1 i 1 i 1 xi 1 x ti 1 Ai 1xi 1 f i 1 Suy ra A x f x x hx O h2 ; i 1 i 1 i 1 i 1 i i 2 Ai 1 xi 1 f i 1 xi 1 xi hxi O h . Thay vào (1.2) ta được h2 h 2 xhx x 2 xxhx xh 2  1 xhx Oh 2  1 x  1 xhxOh 2 . iiiiiii2 2 0 ii 1 iii 2 hay 2 2 1 1 2 3 h xi h 0  1  1 x i O h , tức là 2 2 1 1 1 3 h Ai xi f i h  0 1  2 Ai xi f i O h . Như vậy công thức (3.2) là ổn định và có bậc thấp nhất là bậc một nếu điều kiện sau thỏa mãn: 1 1 1 0  1  2 1. (3.3) Tương tự lược đồ (3.2), xét phương pháp một tựa của nó (one leg variant): 2 xi h A 2 t i  2 x i f  2  t i  . (3.4) Để công thức (3.4) là ổn định và có bậc thấp nhất là bậc một thì các hệ số 2 2 2 0,,  1  2 phải thỏa mãn điều kiện tương tự như điều kiện (3.3), nghĩa là 2 2 2 0  1  2 1. (3.5) Đặt (các giá trị này thay đổi theo từng bước) 2 2 t ti 1  1 2  0 h , A  A t  , f  f t  Với điều kiện (3.5) ta có: 2 2 2 2 2 2 2ti  0112101 t i  t i  t i  t i 2 h  11 t i h  21 t i 2 2 2 2 2 2 2 (0  1  2 )ti 1 (  1 2  0 ) h t i 1  1 2  0 h t . Công thức (3.4) được viết lại thành 53
2 xi h A2  x i  f  hay 2 2 2 2 xi 1 2 x i x i 1 h A  0 x i 1  1 x i  2 x i 1 f  (3.6) 2 2 2 2 2 2 2 EhAx 0 i 1 2 EhAx  1  i EhAx  2  i 1 hf  . Từ (3.2) và (3.6) ta có hệ phương trình như sau: 2 1 2 1 2 1 E h0 Ai 12 E h  1 A i E h  2 A i 1 2  i x x x h , (3.7) 2 2 i 1 2 2 i 2 2 i 1  E h0 A 2 E h  1 A  E h  2 A   trong đó E là ma trận đơn vị cấp n n, i 1  f i ,  f t  . Hệ (3.7) có số phương trình nhiều hơn số ẩn nên ta không thể tìm nghiệm của hệ (3.7) theo nghĩa cổ điển. Để giải hệ (3.7) ta nhân hai vế của nó với ma trận chữ nhật cỡ n 2 n 2 1 2 2 Li 1 E h 0 A i 1 d E h  0 A , với d là hệ số tùy ý. Bằng cách này ta thu được một công thức sai phân hai bước có dạng: Mi 1 x i 1 Px i i Q i 1 x i 1 g i 1 , (3.8) trong đó E h 2  1 A 2E h 2  1 A M L 0 i 1 ; P L 1 i ; i 1 i 1 2 2 i i 1 2 2 E h  0 A 2E h 1 A E h 2  1 A  Q L 2 i 1 ; g L h 2 i . 3. 9 i 1 i 1 2 2 i 1 i 1  E h  2 A  Khi điều kiện (3.3) và (3.5) được thỏa mãn thì công thức (3.2) và (3.6) có cấp chính xác là một, cho nên (3.8) cũng có cấp chính xác là một với mọi giá trị của tham số j i ,i 1,2; j 0,1,2 thoả mãn điều kiện (3.3), (3.5) và với mọi giá trị của d . Điều đó có nghĩa là (3.8) có tất cả năm tham số tự do, chúng quyết định cấp chính xác của công thức. 54
Công thức (3.8) được viết lại như sau 2 2 E h2β 1 A d E h 2 β 2 A x 0i 1 0μ i 1 2 1 2 1 2 2 2 1 E h βA0i 1 2 EhβA 1 i dEhβA 0μ 2 EhβA 2 μ x i E h2βA 1 EhβA 2 1 dEhβA 2 2 EhβA 2 2 x 0i 1 2 i 1 0μ 2 μ i 1 h2 E h 2βA 1 βf 1 βfβf 1 1 dEhβAft 2 2 . 0i 1 0 i 1 1 i 2 i 1 0 μ μ hay 2 1 2 4 1 2 2 2 2 2  1 d E h 2 0 Ai 1 2d 0 A h  0 Ai 1 d  0 A xi 1 2 1 1 2 2 2 2d E h  1 Ai 2 0 Ai 1 d 1 A 2d 0 A x 4 1 1 2 2 2 i h  0  1 Ai 1 Ai d 0  1 A 2 2 1 2 2 1 d E h  1 Ai 1  0 Ai 1 d 2 A d 0 A x 4 1 1 2 2 2 i 1 h  0  2 Ai 1 Ai d 0  2 A 2 1 1 1 h  0 f i 1  1 f i  2 f i 1 df t   4 1 1 1 1 2 h   0 Ai 1  0 f i 1  1 f i  2 f i 1 d 0 A f t   3.8 Khai triển Taylor tại ti 1 ta được: 4h2 8 h 3 fft ft 2 hf 2 hf f fOh 4 i 1 i 1 i 1 i 1 i 12 i 1 6 i 1 2 3 h h 4 fftfthfhfi i i 1 i 1 i 1 f i 1 fOh i 1 2 6 2 2 3 3 B h B h 4 fft  ftBhfBhf i 1 i 1 i 1 f i 1 fOh i 1 2 6 h2 h3 h 4 h 5 y y t y t h y hy y y y 4 y 5 O h6 i i i 1 i 1 i 1 2 i 16 i 1 24 i 1 120 i 1 4h2 8 h 3 AAt At 2 hA 2 hA A AOh 4 i 1 i 1 i 1 i 1 i 12 i 1 6 i 1 2 2 3 3 B h B h 4 AAt  At i 1 BhABhA i 1 i 1 A i 1 AOh i 1 2 6 h2 h 3 AAt At hAhA A AOh 4 i i i 1 i 1 i 12 i 1 6 i 1 Thay vào (3.8’) ta được: 55
2 3 2 2 3 3 2 1 4h 8 h 2 B h B h 1 dEh 20 AhAi 1 2 i 1 A i 1 AdABhA i 1 2  0 i 1 i 1 A i 1 A i 1 2 6 2 6 2 32 2 2 3 3 2 2 2  4 1 4h 8 h 2 B h B h h 0 AhAi 1 2 i 1 A i 1 Ad i 1  0 ABhA i 1 i 1 A i 1 A i 1  2 6 2 6  2 3 4 5 h h h 4 h 5 xi 1 hx i 1 x i 1 x i 1 x i 1 x i 1 2 6 24 120 2 3 2 3 2 1 4h 8 h 1 4 h 8 h 2 2dEh 1 AhAi 1 2 i 1 A i 1 A i 1 2  0 AhA i 1 2 i 1 A i 1 A i 1 2 6 2 6 2 2 3 3 2 2 B h B h d1 2 d  0 Ai 1 BhA i 1 A i 1 A i 1 2 6 2 3 2 3 4 1 1 4h 8 h h h h 0  1 AhAi 1 2 i 1 A i 1 AAhA i 1 i 1 i 1 A i 1 A i 1 2 6 2 6 2 2 3 3 2  2 2 B h B h d0  1 Ai 1 BhA i 1 A i 1 A i 1  2 6  2 3 4 5 4h 8 h 16 h 4 32 h 5 xi 1 2 hx i 1 x i 1 x i 1 x i 1 x i 1 2 6 24 120 2 3 2 2 1 4h 8 h 1 d E h 1 Ai 1  0 A i 1 2 hA i 1 A i 1 A i 1 2 6 B2 h 2 B 3 h 3 d2 d  2 A BhA A A 2 0 i 1 i 1 i 1 i 1 2 6 2 3 2 3 4 1 1 4h 8 h h h h 0  2 Ai 1 2 hA i 1 A i 1 A i 1 A i 1 hA i 1 A i 1 Ai 1 2 6 2 6 2 2 3 3 2 2 2 B h B h d0  2 Ai 1 BhA i 1 A i 1 A i 1 x i 1 2 6 2 3 2 3 2 1 h h 1 h h h 0 fhfi 1 i 1 f i 1 f i 1  1 fhf i 1 i 1 f i 1 f i 1 2 6 2 6 2 2 3 3 1 B h B h  2fi 1 d f i 1 Bhf i 1 f i 1 f i 1  2 6  56
2 3 2 3 4 1 h h 1 h h h 0 AhAi 1 i 1 A i 1 A i 1  0 fhf i 1 i 1 f i 1 f i 1 2 6 2 6 2 3 1 h h 1 1 fi 1 hf i 1 f i 1 f i 1  2 f i 1 2 6 2 2 3 3 2 2 3 3 2 B h B h B h B h  d0 Ai 1 BhA i 1 A i 1 A i 1 f i 1 Bhfi 1 f i 1 f i 1 . 2 6 2 6  Hệ số của h0 1 d xi 1 2 2 d x i 1 1 d x i 1 0 đúng với mọi d . Hệ số của h1 : 1 d 2 xi 1 2 2 d x i 1 0 đúng với mọi d . Hệ số của h2 : 1 1 dx 2 21 2 dAx  2 2 2 dx .  1 2  1 ddAx  2 2  2 i 1 0011 i i 2 i 1101011 i i 1 1 2 2 1 1 1 2  0 d  2 d  0 Ai 1 x i 1  0  1  2 d f i 1 1d x 1  1  1 d  2  2  2 A x  1  1  1 d f . i 1 012 01211012 i i i 1 Nếu điều kiện (3.3) và (3.5) được thoả mãn thì hệ số của h 2 bằng 0, suy ra công thức (3.8) có cấp chính xác bằng một. Hệ số của h3 : 8 1 dx 21 AdAx 2  2 2 2  1 .2 AdBAx 2  2 6 i 1 01 i 011 i i 01 i 011 i i 1 2 2dx 1 2  1 ddAx  2 2  2  1 2  1 .2 ddBAx  2 2  2 6 i 1101 01110 i i 1 011 i i 1 2 2 1 1 20 dB  2 dB  0 Ai 1 x i 1 2  0  1 dB f i 1 1 2 1 2 1 2 1 dx i 1 2 0 2 dAx  01111111 i i  dAx  i i  dBAx  111 i i 1 2 2 1 1 20 dB  2 dB  0 Ai 1 x i 1 2  0  1 dB f i 1 1 2 1 2 1 d xi 1 2 0 2 d  0  1 d  1 A i 1 x i 1 1 2 1 2 2 1 1 1 dB  1 2  0 dB  2 dB  0 Ai 1 x i 1 2  0  1 dB f i 1 1 1 1 d xi 1 2 0  1 dB A i 1 x i 1 A i 1 x i 1 f i 1 57
Từ (3.1) suy ra x t A t x t A t x t f t * nên hệ số của h3 được viết lại như sau: 1 1 1 d xi 1 2 0  1 dB x i 1 . Vậy hệ số của h3 bằng 0 khi và chỉ khi 1 1 2 2 1 d 20  1 dB ; B  1 2  0 . (3.10) Hệ số của h4 : 16 1 d x 4 21 2 d  2 A .2 x 4  1 2 dB  2 A .2 x 24 i 1 0 0 i 1 i 1 0 0 i 1 i 1 2 2 2 41 dB 2  2 A x  1 d  2 A 2 x 1 d x 4 0 0 i 1 i 1 0 0 i 1 i 124 i 1 1 1 2  1 d  2 2 d  2 A . x  1 4  1 dB  1 2 dB  2 A x 101 01110 i 2 i 2 011 i i 11 1 1 2 2 2 2 1 1 2 2 2 1 4  0 dB 1 dB  0 Ai 1 x i 1  0  1 d  0  1 A i 1 x i 1 2 2 11 2 2 1 2 2 1 1 2 2 2 20 dB  2 dB  011 Ai x i  020211  d   A i x i 2 2 11 1 1 2 1 2 20  1 dB fi 1  0 d  0 A i 1 f i 1 . 2 2 Từ (*) suy ra x 4 t Atxt Atxt Atxt Atxt ft Atxt 2 Atxt AtAtxt ft ft Atxt 2 Atxt Atxt2 Atft ft Vậy hệ số của h4 có thể viết lại như sau 58
14 14 14 14 d 1 A x d 1 A x d 1 A2 x d 1 A f 24 i 1 i 1 24 i 1 i 1 24 i 1 i 1 24 i 1 i 1 14 d 1 f 4 1 4d 2 A2 x 4 1 4d 2 A f 24 i 1 0 0 i 1 i 1 0 0 i 1 i 1 1 2 1 2 2 1 2 2 2 2 8 0 4dB 0 Ai 1 xi 1 4 0 dB  0 Ai 1 xi 1  0 d  0 Ai 1 xi 1 1 1 1 1 2 2 2 1 1 1 1 2 2 1  0 d1 d 0 Ai 1 xi 1 1  0 d1 d 0 Ai 1 fi 1 2 2 2 2 1 1 1 2 1 1 1 1 2 2 2 1 4 0 dB 2 2dB 0 Ai 1 xi 1 1 2 0 dB dB  0 Ai 1 xi 1 2 2 1 1 2 2 1 1 2 2 2 1 1 1 1 2  0 1 d 0 1  0 1 d 0 1 Ai 1 xi 1 2 0 1 dB f i 1 2 2 71 2 1 1 1 1 2 2 1 d 40 dB  0 (  1 2  0 dB dB  0 ) Ai 1 x i 1 12 2 2 71 2 1 1 2 2 7 1 2 1d 80 4 dB  0 (  1 4  0 dBdBAx  1 2  0 )i 1 i 1 1 d 4  0 4 d  0 6 12 2 2 1 2 1 11 1 221122112 2 2 0 d  0  10 d  100101020 d  d d 2 Ai 1 x i 1 2 2 71 2 1 1 1 1 2 2 1 2 1d 40 4 d  0  1  0 d  1 d  0  0 d  0 Ai 1 f i 1 12 2 2 7 11 1 1 2 1 d 1 2  0 dB fi 1 12 2 2 7 1 1 2 1 d 40  1 dB Ai 1 x i 1 A i 1 x i 1 A i 1 f i 1 f i 1 5 7 1 21 1 1 2 1 d 20 2 d  0  1 Ai 1 x i 1 . 12 2 2d 2 1 7 1 1 2 1 d 40  1 dB ; 4 6 Hệ số của h bằng 0 khi 7 1 1 1 d 21 2 d  2  1 d  2 . 120 0 2 1 2 1 1 1 Từ điều kiện (1.10) ta có 20 1 d 1 dB . Thay vào hệ trên ta được 5 1 d  1 dB2 2dB 6 1 3.11 5 1 d  1 d 2 6 1 1 59
Hệ số của h5 : 32 4 1 d x 5 21 d  2 A . x 2 2  1 dB  2 A .2 x 120i 1 0 0 i 1 3 i 1 0 0 i 1 i 1 2 2 1 41 dBAx 2  2 .2  1 d  2 Ax 2 .2 8  1 dBAx 3  2 0 0110 i i 011 i i 3 0 011 i i 2 2 1 1 2 21 dB  2 A A x 2 1 d x 5  1 2  1 d  2 2 d  2 A . x 0 0 i 1 i 1 i 1120 i 1 1 0 1 0 i 16 i 1 1 1 1 4  1 dB  2 2 dBAx  2 .  1 8  1 dB 2  2 2 dBAx 2  2 101 01110 i 2 i 2 1 011 i i 1 1  1 d  2  2 A 2 x  1 16  1 6 d  2 2 dB 3  2 A x 010111 i i 6 1 0 1 011 i i 1 31  1 2dB  2  2 A A x 8  1 dB 3  2 dB 3  2 A x 0 1 0 1 i 1 i 1 i 16 0 2 0 i 1 i 1 1 31  1 2dB  2  2 A A x 8  1  1 dB 3 f 0 2 0 2 i 1 i 1 i 16 0 1 i 1 2 2 21  1  1A f 2  1  1  1 A f dB  2 A f A f . 0 0111 i i 0 0111 i i 011 i i i 11 i Từ (*) và ( ) suy ra x 5 t Atxt Atxt 2 Atxt 2 Atxt Atxt f t Atxt 3 Atxt 3 Atxt Atxt f t Atxt 3 Atxt 3 AtAtxt ft AtAtxt Atxt ft f t Atxt 3 Atxt 4 AtAtxt Atxt2 3 Atft Atft f t . Hệ số của h5 được viết lại như sau: 1 3 4 1 3 dAx 1 dAx 1 dAAx 1 dAx 1 2 dAf 1 4ii 11 4 ii 11 4 iii 111 4 ii 11 4 ii 11 1 1 8 dAf 1 df 1 1 dAAxAxf  2 4i 11 i 4 i 1 3 00111111 i i i i i i 2 2 4 21 dB  2 A A x f 2 4  1 dB 2  2 A x 2 1 d  2 A 2 x 0 0 i 1 i 1 i 1 i 1 0 0 i 1i 1 0 0 i 1 i 1 1 2 2 81 dB 3  2 A x 2 2  1 dB  2 A A x 3 0 0 i 1 i 1 0 0 i 1 i 1 i 1 1 1 2  1 d  2 2 d  2 A A x A x f 6 1 0 1 0 i 1 i 1 i 1 i 1 i 1 i 1 1 1 1 4  1 dB  2 2 dBAAxf  2  1 8  1 dB 2  2 2 dB 2  2 Ax 2 101 01111 i i i i 2 10 1 011 i i 60
1 1  1 d  2  2 A 2 x  1 16  1 6 d  2 2 dB 3  2 A x 010111 i i 6 1 0 1 011 i i 1 31  1 2dB  2  2 AAx 8  1 dB 3  2 dB 3  2 Ax 01 01111 i i i 6 0 2 011 i i 1 31  1 2dB  2  2 A A x 8  1  1 dB 3 f 0 2 0 2 i 1 i 1 i 16 0 1 i 1 2 2 21  1  1 dB  2 A f 2  1 2  1  1 dB  2 A f 0 01 011 i i 0 01 011 i i 1 1132 1 1 12 32 1 13232 d1 800101 dB   16  ddB  2  0 8  02011 dBdBAx   i i 4 3 6 6 31 2 2 1 1 1 2 2 2 2 d 1 2 40 dB  0  1 8  0 dB  1 2 dB  0 Ai 1 x i 1 4 2 8 1 2 1 2 12 2 2 d 1 0 d  0 4 2  0 dB  0 2 2  0 dB  0 3 1 1 1 4  1 2dB 2 dB  2  1 2  1 d  2 2 d  2 2 1 0 1 0 6 1 0 1 0 31  1 2dB  2  2 3  1  1 2 dB  2  2 A A x 01 01 02 02 i 111 i i 1 81 2 12 2 2 1 1 1 2 2 d 1 0 d  0 2  0 d  0  1 2  0 d  1 2 d  0 4 3 6 11222 3 1211 1 1 22 010111 d   Ai x i d 1 4 2  00 dB   10 2  dB  10 dB  4 2 2 2 1 1 1 2 2 0  0  1 dB  0 Ai 1 f i 1 2 1 81 2 1 1 1 2 2 1 1 1 2 d 1 0 d  0  1 2  0 d  1 2 d 0 2  0  0  1dB  0 Ai 1 f i 1 4 3 6 1 1 1 1 3 d 1 80  1 dB fi 1 . 4 6 1 81 1 1 1 3 d 1 0  1 dB Ai 1 x i 1 A i 1 x i 1 f i 1 4 6 6 6 2 1 71 1 1 7 2 1 1 1 1 2 2 2 d1 01  d  0 2  001   d  1 dB  0 Ai 1111 x i A i f i 4 3 6 3 6 2 161 2 1 7 2 1 1 2 1 2 2 d 1 0  1 d  0  0 d  1 d  1 dB  0 Ai 1 A i 1 x i 1 3 3 3 2 6 2 3 1 1 1 7 2 1 1 1 1 2 2 d 1 60  1 d  0 2  0 2  0  1 dB  1 2 dB  0 Ai 1 f i 1 4 2 3 2 61
Hệ số của h5 bằng 0 khi và chỉ khi 1 8 1 1 d 1 1  1 dB 3 ; 4 60 6 1 6 1 7 1 72 1 d 1 1  1 d  2 2  1  1  1 d  2 dB  2 ; 4 30 6 1 3 0 0 0 1 6 1 0 16 2 72 1 1 d 1 1  1 d  2  1 d  2 d  2 dB  2 ; 30 3 1 3 0 0 2 1 6 1 0 3 1 72 1 1 1 2 1 1 1 2 2 d 1 60  1 d  0 2  0 2  0  1 dB  1 2 dB  0 . 4 2 3 2 Từ (3.10) và (3.11) ta có 5  1 d 1 dB 2 2 dB ; 1 6 1 1 1 1 1 1 1 1 (d 1)  1 dB ( d 1) d 1 dBdBdB 2 02 2 1 2 2 12 2 2 5 1 d 1 dB2 dB ; 12 2 2 2 2 3 2 d1 dB 2 dB dB  1 dB 2 dB ; 2 2 2 2 2 3 2 20 B  1 2 dB  0 dBdB  1 dBdB . Suy ra hệ số của h5 bằng 0 khi và chỉ khi Bd 2 3B B2 0 3.12 1 5 2 1 3 d d 1 B B B 0 12 12 2 Nếu các tham số thoả mãn điều kiện (3.3), (3.5), (3.10) và (3.11) thì sai số địa phương của công thức (3.8) là O( h5 ), lúc này các hệ số là nghiệm của hệ phương trình: 1 1 1 0  1  2 1; 2 2 2 0  1  2 1; 1 1 d 1 1 2  0 dB ; 5 d 1 1 d  2 ; 6 1 1 5 d 1  1 dB 2 2 dB . 6 1 Nghiệm của hệ này là 62
Hoặc 5 1 1 d 0,1 ,  1 ,  1 (3.13) 16 0 12 2 12 2 2 2 với các tham số 1,,,  0  2 B tuú ý; Hoặc 5  1 d 1 dB 2 2 dB ; 1 6 2 2 1 BB 2 ; 1 1 1  1 d 1 dB 2 dB ; 0 12 2 2 11 1 1 3 (3.14)  1 d dB 2 dB; 2 12 12 2 2 1 1  2 B 2 ; 0 2 2B 1 3  2 BB 2 1. 2 2 2 với tham số B tuỳ ý, d ≠ 0. Nhận xét 3.1 Khi các tham số thoả mãn điều kiện (3.13) thì công thức (3.8) trở thành: 25 2 1 2 2 1 5 1 RxREhAxRiii 1 1 1(2 iii ) 1 E hAxhR iiiiii 1 1 1 f 1 f f 1 (3.15) 6 12 12 6 12 1 trong đó R E h2 A . i 112 i 1 Sau khi chia hai vế của công thức (3.15) cho Ri 1 thì ta thu được công thức nổi tiếng Numerov có sai số địa phương là O( h6 ). Nhớ lại rằng t ti 1 hB . Đặt 2 1 1 2 1 S i 1 E h d 1 dB dB Ai 1 ; 12 2 2 2 1 2 1 U i 1 E h B B A 2 2 Điều kiện (1.14) vẫn còn hai tham số tự do và ta có phương trình sai phân sau: 63
2 2 Si 1 dU i 1 y i 1 2 5 2 2 2 Si 12 Eh d 1 dB 2 dBAdU i i 1 2 EhB 2 BAy  i 6 (3.16) 21 2 2 1 2 3 SEhi 1( 1 11 d 6 dB 18 dBAdU i 1 i 1 EhB B 1 Ay i 1 12 2 2 2 h Si 1 1 f i dU i 1 f , trong đó hệ số của toán tử 1  fi  được xác định bởi công thức (3.14). Sai số địa phương của công thức (3.16) là O( h5 ). Bởi vì chúng ta vẫn còn hai tham số tự do là B và d , nên vẫn còn có ý nghĩa xét phương trình (3.12). Nghiệm của hệ mở rộng là: Hoặc 5 1 1 d 0,1 ,  1 ,  1 và 2,,,  2  2 B tuỳ ý. 16 0 12 2 2 1 0 2 Với nghiệm này ta lại thu được công thức Numerov. Hoặc 1 1 1 2 2 2 d 1, B 0,1 0,  0 0,  2 1,  1 0,  0 0,  2 1. Hoặc 1 1 1 2 2 2 d 1, B 2,1 0,  0 1,  2 0,  1 0,  0 1,  2 0. Hai trường hợp này đều cho Mi 1 P i Q i 1 g i 1  0 nên không có ý nghĩa. Hoặc 5 1 1 1 d 1, B 1,1 d ,  1 d 1 ,  1 d 1 , 16 6 0 12 2 12 (3.17) 2 2 2 1 1,  0 0,  2 0. Khi các hệ số thoả mãn điều kiện (3.17) thì công thức (3.8) có sai số địa phương là O( h6 ). Không thể chọn được các giá trị tương ứng của d để tăng độ chính xác cao hơn được nữa. 3.1.3. Tính chất ổn định Với phương trình thử x k2 x ta có 64
2 xi 1 2 R mm  x i x i 1 0 , trong đó  hk , hàm phân thức Rmm được gọi là hàm ổn định của phương pháp. 2 Khoảng hội tụ của phương pháp là 0,0 nếu 2 2 2 Rmm  1  0  0 . 2 Phương pháp là ổn định - P nếu Rmm  1 với mọi số thực . Lưu ý rằng phương pháp Numerov có bậc hội tụ cao hơn mọi lược đồ sai phân hai bước, tuy nhiên khoảng hội tụ khá hẹp, chỉ là 0,6 . Trong mục này ta sẽ nghiên cứu các tính chất ổn định của các phương pháp (3.8) với các hệ số (3.17) đem lại lược đồ phụ thuộc vào một tham số d 0 có bậc bốn h2 1 V E d 1 A ;   f d 1 f 10 2 d f d 1 f . i 112 i 1 i12 i 1 i i 1 Ta có 2 2 t ti 1  1 2  0 h t i 1 h t i h2 d 1 AAffL ,, E AdEVdE i  i i 1 i 1 i 1 12 Và với điều kiện (3.17) thì (3.8) được viết lại thành 2 2 2 h 2 h VdExi 1 i 1 V i 1 2 E 10 2 dAdEhAxV i 2 i i i 1 E d 1 A i 1 dEx i 1 12 12 2 h Vi 1 f i df i . 1.18 Áp dụng công thức (3.18) cho phương trình thử y k2 y ta được 2 2d 1  2 d 1 2  2 E dE xi 1 E 2 E 10 2 d d 2 E  x i 12 12 12  2 d 1  2 E E d 1 dE x 0. i 1 12 12 Hay 65
2 2  d 1  10 2 d 2 E 2 E d 2 E  1 12 12 xi 1 2 2 x i x i 1 0 2  2 d 1 E dE 12 với  hk . Như vậy công thức (3.18) có hàm ổn định tương ứng là 2 2  d 1  10 2d 2 E 2E d 2E  12 12 2 1 R 22 R 22  2 2  2 d 1 E dE 12  2  4 1 5 d 3 144 3.19  2  4 1 d 1 6 144 Ta có 2 RR22  1   khi và chỉ khi 4 2d 0 4  2 9 4 2d 2 0  R 1 8 4 2d d . 144 6 0 4 36 144 9 9 Như vậy những công thức ứng với d là ổn định - P. Với d các phương 4 4 pháp có khoảng hội tụ được mở rộng hơn so với công thức Numerov. Bảng dưới đây đưa ra một số ví dụ để so sánh. Nhận xét rằng với d 0 ta nhận được phương pháp Numerov. d  2 d  2 d  2 d  2 d  2 0 6 0,5 4 4 7 6,58 1 6 6 5 7, 41 1,5 12 12 3 8,78 2 12 2 Một số đồ thị của R22 () với d 2,5 được chỉ ra trong Hình 3.1. Đồ thị phù hợp với lý thuyết đã trình bày ở trên. 66
Hình 3.1 3.1.4. Thử nghiệm số Trong phần này ta tính toán theo lược đồ 3.18 với các giá trị khác nhau của d để được các kết quả số trên hai ví dụ kiểm tra. Chúng ta giả thiết rằng khi thực hiện các lược đồ sai phân ta đã biết các giá trị chính xác của dữ liệu đầu vào x0 và x1 . Ví dụ 3.1 Xét phương trình 2a a2 x 3 4 x, t t 67
a a 2 a với nghiệm tổng quát là x t C et C e t e t dt, trong đó CC, là hai hằng số 1 2 1 2 tuỳ ý. Cho C1 1, C 2 0, t  1, 10 , a 20. Xây dựng trên đoạn [1,10] một lưới với bước 0.2 9 lưới là h , j 1,2,3 và cho sai số là er e 2 y, N . j j j Nj j 2 hj Mã nguồn và Chương trình được cho trong Phụ lục. Kết quả tính toán được cho trong Bảng 1 dưới đây. d h = 0.1 h = 0.05 h = 0.025 2 6.3.10 3 5.8.10 4 4.5.10 5 0 2.9.10 4 3.1.10 5 2.4.10 6 3 4 5 2 5.2.10 5.2.10 4.10 3 7.7.10 3 7.9.10 4 6.2.10 5 5 1.2.10 2 1.3.10 3 1.10 4 Bảng 1. Sai số cho x của ví dụ 1 với các giá trị khác nhau của h và d . Ví dụ 3.2 Xét hệ phương trình x k2 x 2 ak cos kt , x 0 1, x 0 0; y kyak2 2sin kty , 0 0, y 0 kat ,  0,40,  với nghiệm chính xác là x cos kt at sin kt , y sin kt at cos kt . Trường hợp k 1 và a 0,0005 đã được xét trước đây. Trong Ví dụ này, ta xét giá trị k 5 . Ta đưa vào đại lượng  t x2 t y 2 t 1 at 2 . Xây dựng trên đoạn 0, 40  một lưới với bước 68
hj , j 1, 2,3,4 er j  40  h N j , 4 j j j trong đó Nj 40 4 h j . Ta có công thức sai phân tương tự  i x2 y 2 , trong đó các giá trị x, y được hj i i i i tính từ công thức (3.18). Với các dữ liệu đầu vào thỏa mãn điều kiện (3.17), phương pháp (3.18) cho bước tích phân (với một số giá trị của tham số d ) rộng hơn đáng kể. Ta cũng lưu ý rằng khi d 1 lược đồ sai phân là ổn định với  2 0;6 5 1 và khi d 2 lược đồ sai 9 phân là ổn định với  2 0;12 và là P-ổn định với d . 4 Mã nguồn và Chương trình được cho trong phần Phụ lục. Kết quả tính toán cho các giá trị của d và các bước hj khi k 5, a 0.0005 được cho trong Bảng 2. d j 1 j 2 j 3 j 4 1 * 1.9.10 2 9.8.10 5 3.8.10 7 2 * 4.2.10 2 3.4.10 4 1.3.10 6 3 3 2 4 6 5.2.10 2.3.10 5.9.10 2.3.10 4 2 5.6.10 2.8.10 3 8.3.10 3 3.2.10 6 5 1 1.0.10 1.3.10 2 1.0.10 3 4.2.10 6 Bảng 2. Chú thích Dấu * là kí hiệu sai số lớn hơn 106 . 3.2. Phương pháp không cổ điển giải số hệ phương trình vi phân phi tuyến cấp hai Các kết quả trong 3.1 có thể mở rộng cho hệ phi tuyến bậc hai như sau. Xét hệ phương trình phi tuyến bậc hai xtgxttx () ((),), (0) xx0 , (0) xt 0 ,  0,1 . (3.20) Tương tự lược đồ (3.2), ta viết lược đồ sai phân sau đây cho bài toán (3.20): 2 1 1 1 xi 1 2 xxh i i 1  0 gxt i 1 , i 1  1 gxt i , i  2 gxt i 1 , i 1 . (3.21) 69
Để giải hệ (3.21) ta áp dụng phương pháp Newton để nhận được công thức j Mij x i 1 l ij với 2 1 j Mij E h0 J x i 1, t i 1 và l 2 yyh 2 1 gxtj 1 , Jxtx j 1 , j 1  1 gxt j 1 ,  1 gxt , . ij i i 1 0 i i 1 i 1 i 1 i 1 1 i i 1 2 i 1 i 1 Ở đây chỉ số trên j kí hiệu số bước lặp và J(,)(,) x t  x g x t . Viết lại hệ trên và thao tác tương tự như (3.6), ta nhận được hệ phương trình đại số tuyến tính M lij ij x j , (3.22) i 1 E mij trong đó 2 mij 2 x i x i 1 h g ( x i , t i ) . Hệ (3.22) không có nghiệm cổ điển. Ta cố gắng tìm “nghiệm” của hệ (3.22) tương tự như trong 3.1, tức là nhân hai vế của (3.22) với ma trận Mij dE . Kết quả ta nhận được hệ không suy biến 2 j Mij dE x i 1 M ij l ij dm ij . (3.23) 0 Ở đây ta chọn xi 1 x i làm xấp xỉ ban đầu. Mặc dù các kết quả lý thuyết cho hệ phi tuyến chưa thật sự trọn vẹn, các ví dụ trong [4] chỉ ra rằng, lược đồ (3.23) tỏ ra khá hiệu quả. 70
KẾT LUẬN Luận văn đã trình bày một phương pháp giải số không cổ điển giải bài toán Cauchy cho hệ phương trình vi phân bậc nhất do M. V. Bulatov mới đề nghị gần đây (2003-2008) và một phương pháp giải số bài toán Cauchy cho hệ phương trình vi phân bậc hai không cổ điển của M. V. Bulatov và G. V. Berghe (2009). Mặc dù các kết quả lý thuyết cho hệ phi tuyến bậc hai cũng như các phương pháp tổng quát không đối xứng cho cả hệ tuyến tính và hệ phi tuyến chưa thật sự trọn vẹn, các ví dụ trong [4] chỉ ra rằng, lược đồ (3.23) tỏ ra khá hiệu quả. Vì vậy phương pháp do Bulatov (và Berghe) đề xuất đáng đựoc quan tâm và tiếp tục phát triển. Trong quá trình trình bày các kết quả trên, chúng tôi đã phải tính toán công phu, tỷ mỷ và cẩn thận trong chứng minh lý thuyết (phân tích các hàm nhiều biến vào chuỗi Taylor, diễn giải chi tiết các công thức), đồng thời chúng tôi cũng lập trình trên MATLAB và thực hiện tính toán trên máy nhằm kiểm nghiệm lại bốn ví dụ trong [4] và [11]. Tiếc rằng không có đủ thời gian để chúng tôi có thể sử dụng phần mềm tính toán Maple để tính tự động đạo hàm cũng như sử dụng cơ sở Grobner để giải các hệ phương trình đa thức nảy sinh trong các lược đồ mở rộng các kết quả của [4] và [9]- [11]. Hy vọng những vấn đề này sẽ được tiếp tục nghiên cáu trong thời gian tới. 71
TÀI LIỆU THAM KHẢO [1] Phạm Kỳ Anh (1996), Giải tích số, Nhà xuất bản Đại học Quốc gia Hà Nội. [2] Tạ Duy Phượng (2009), Một số Chương của Giải tích số và thực hành tính toán, Giáo trình Cao học, Viện Toán học, Hà Nội. [3] U. M. Ascher, L. R. Petrold (1998), Computer Methods for Ordinary Differential Equations and Differential-Algebraic Equations, SIAM. [4] M. V. Bulatov, G. V. Berghe (2009), Two-step fourth order methods for linear ODEs of the second order, Numer. Algorithms, 51, No4, pp. 449-460. [5] E. Hairer, S. P. Nørsett, G. Wanner (1993), Solving Ordinary Differential Equations I Nonstiff Problems, Springer-Verlag. [6] J. D. Lambert (1991), Numerical Methods for Ordinary Differential Systems The Initial value Problem, John Wiley & Son Ltd. England. [7] D. Quiney (1987), An Introduction to the Numerical Solution of Differential Equations, John Wiley & Son Inc., England. [8] J. Stoer, R. Burlisch (2002), Introduction to the Numerical Analysis, Springer- Verlag, New York. [9]. M. В. Булатов (2003), О построении 1-стадийного L-устойчивого метопа второго порядка, Дифференциальные уравнения, том 39, No4, с.554-556. [10]. M. В. Булатов (2005) Построение неклассичиских многошаговых схем для линейнынх оду , Доклады академии наук, том 404, No1, с.11-13. [11]. M. В. Булатов (2008), О построении неклассических разностных схем для обыкновенных дифференциальных уравнений, Дифференциальные уравнения, том 44, No4, с.546-557. 72
MỤC LỤC Lời nói đầu 1 1. Kiến thức chuẩn bị 4 1.1. Bài toán Cauchy giải hệ phương trình vi phân 4 1.2. Giải số bài toán Cauchy 5 1.2.1. Quy tắc cầu phương cơ bản và giải số phương trình vi phân 5 1.2.2. Phương pháp Runge-Kutta 10 1.2.3. Phương pháp cổ điển đa bước 13 1.3. Mô hình thử và ổn định của phương pháp số 14 1.3.1. Mô hình thử 14 1.3.2. Ổn định của phương pháp Euler 15 1.3.3. Ổn định của phương pháp Runge-Kutta 17 1.3.4. Ổn định của phương pháp đa bước 19 1.3.5. Ổn định của phương pháp sai phân hữu hạn 19 2.Về một phương pháp không cổ điển giải số hệ phương trình vi phân cấp một 21 2.1. Phương pháp không cổ điển giải số hệ phương trình vi phân phi tuyến cấp một 21 2.1.1. Phương pháp tổng quát 21 2.1.2. Phương trình thử 25 2.1.3. Trường hợp đặc biệt 26 2.1.4. Thử nghiệm số 27 2.2. Phương pháp không cổ điển giải số hệ phương trình vi phân tuyến tính cấp một 28 2.2.1. Phương pháp một bước 28 2.2.2. Phương pháp đa bước 35 i 73
3. Về một phương pháp không cổ điển giải số hệ phương trình vi phân cấp hai 51 3.1. Phương pháp không cổ điển giải số cấp hai 51 3.1.1. Phương pháp cổ điển 51 3.1.2. Lược đồ sai phân mới 52 3.1.3. Tính chất ổn định 64 3.1.4. Thử nghiệm số 67 3.2. Phương pháp không cổ điển giải số hệ phương trình phi tuyến cấp hai 69 Kết luận 71 Tài liệu tham khảo 72 Phụ lục 72 ii74