Tính vững, ổn định và hội tụ của phương pháp sai phân hữu hạn cho phương trình nhiệt

Nội dung text: Tính vững, ổn định và hội tụ của phương pháp sai phân hữu hạn cho phương trình nhiệt

1. Tạp chí Khoa học và Công nghệ, Số 43B, 2020 TÍNH VỮNG, ỔN ĐỊNH VÀ HỘI TỤ CỦA PHƯƠNG PHÁP SAI PHÂN HỮU HẠN CHO PHƯƠNG TRÌNH NHIỆT NGÔ NGỌC HƢNG Trường Đại học Công nghiệp Thành phố Hồ Chí Minh ngongochung@iuh.edu.vn Tóm tắt. Trong bài báo này tôi sẽ bàn về phƣơng pháp số để giải phƣơng trình nhiệt với điều kiện ban đầu và điều kiện biên Dirichlet. Xấp xỉ của các đạo hàm bằng phƣơng pháp sai phân hữu hạn giữ vai trò quan trọng trong phƣơng pháp số trong lĩnh vực phƣơng trình đạo hàm riêng, đặc biệt các bài toán biên. Việc nghiên cứu tính nhất quán và tính ổn định của nghiệm xấp xỉ là cần thiết. Vì có các tính chất này, nghiệm xấp xỉ mới đảm bảo hội tụ về nghiệm chính xác. Ví dụ số cũng sẽ đƣợc thực hiện để minh họa cho các kết quả lý thuyết. Từ khóa. Phƣơng trình nhiệt, phƣơng pháp sai phân hữu hạn, tính vững, tính ổn định. CONSISTENCY, STABILITY AND CONVERGENCE OF FINITE DIFFERENCE SCHEMES ON THE HEAT EQUATION Abstract. This paper deal with a numerical method for the solution of the heat equation together with initial condition and Dirichlet boundary conditions. The approximation of derivatives by finite differences plays a central role in finite difference methods for the numerical solution of differential equations, especially boundary value problems. The consistency and the stability of the schemes are described. Futhermore, numerical simulations are performed to illustrate the accuracy and stability of the regularized solution. Keywords. Heat equation, finite difference method, consistency, stability. 1 GIỚI THIỆU Trong thực tế, nhiều vấn đề trong vật lý nhƣ phƣơng trình truyền nhiệt, phƣơng trình sóng, phƣơng trình Poisson và phƣơng trình Laplace đƣợc mô hình hóa bằng các phƣơng trình đạo hàm riêng. Một số phƣơng trình đạo hàm riêng này có nghiệm chính xác trong những miền đặc biệt. Nhƣng nói chung việc xác định nghiệm chính xác của các phƣơng trình đạo hàm riêng trên miền bất kỳ gặp nhiều khó khăn. Do đó, việc nghiên cứu phƣơng pháp tính số để tìm nghiệm gần đúng là rất quan trọng. Phƣơng pháp sai phân hữu hạn là một trong những phƣơng pháp dùng để tìm nghiệm xấp xỉ của phƣơng trình vi phân. Bằng việc phân hoạch miền xác định thành hữu hạn các lƣới nhỏ. Nghiệm xấp xỉ đƣợc tính tại các điểm lƣới của miền xác định [1]. Bài viết này liên quan đến phƣơng pháp sai phân hữu hạn cho phƣơng trình nhiệt trong thanh vật liệu có chiều dài L , phƣơng trình có dạng nhƣ sau ut  u xx fxt , , xt , (0,)(0,], L T (1) với điều kiện ban đầu u( x ,0) g (x ), x [0,L ], (2) và điều kiện biên Dirichlet u(0,tht ) ( ), t [0, T ], (3) uLt(,) kt (),t [0,],T (4) trong đó hàm chƣa biết u( x , t ) là giá trị nhiệt độ tại vị trí ở thời điểm t, ut ( x , t )và uxx ( x , t )tƣơng ứng là đạo hàm riêng cấp một và cấp hai của hàm u( x , t ) theo biến thời gian t và không gian x . Hằng số  là © 2020 Trƣờng Đại học Công nghiệp Thành phố Hồ Chí Minh
2. 106 TÍNH VỮNG, ỔN ĐỊNH VÀ HỘI TỤ CỦA PHƢƠNG PHÁP SAI PHÂN HỮU HẠN CHO PHƢƠNG TRÌNH NHIỆT độ dẫn nhiệt của vật liệu. Hàm số g( x ) là phân bố nhiệt độ tại thời điểm ban đầu. T là thời điểm cuối. Các hàm số theo biến thời gian h( t ) và k( t ) mô tả dòng nhiệt tại hai biên. Ở đây tôi giả định là vật liệu đồng nhất và bề mặt của thanh đƣợc cách nhiệt để nhiệt chỉ truyền theo hƣớng x. Giả sử thêm là bài toán trên là chỉnh và có nghiệm duy nhất u( x , t ) . 2 PHƯƠNG PHÁP SAI PHÂN HỮU HẠN Đầu tiên, ta tiến hành rời rạc hóa miền liên tục [0,L ] [0, T ] thành một tập Nx N t điểm lƣới. Ở đây tôi sử dụng các điểm chia cách đều nhau, nghĩa là ta chia miền của x thành tập các điểm cách đều nhau L xi i x,, x ,,, iN 0  x (5) Nx tƣơng tự, đối với miền của t đƣợc chia nhƣ sau T tjjj i ,. ,,, jN 0  t (6) Nt j Ta ký hiệu ui là giá trị của u( x , t ) tại điểm lƣới (xi , t j ). Ta xấp xỉ đạo hàm riêng ut( x i , t j ) của phƣơng trình (1) nhƣ sau uj 1 u j uxt( , ) i i ,j  0, , N 1, (7) tij t t trong đó chú ý, từ điều kiện ban đầu (2) ta có 0 u(xugi ,0) i ( xi i ),, 0,  N x . (8) Tƣơng tự, đối với đạo hàm riêng cấp hai uxx( x i , t j ) ta có xấp xỉ nhƣ sau uj 2 uu jj u (x,t ) i 1 ii 1 , i  1, , N 1, j 0,  ,N 1, (9) xx i j 2x x t và theo điều kiện Dirichlet, ta có j u(0,tuhtj ) 0 ( j ), j 0,  , N t 1, (10) và uLtu( , ) j kt (),j 0, , N 1. (11) jNjx t Vậy ta có công thức xấp xỉ cho (1) nhƣ sau uuj 1 j u j 2 uu jj i i  i 1 ii 1 j 2 fi t x ,i 1,,Nx 1, j 0,  , N t 1 (12) uj ht(), uj k ( t ), u0 g( x ) 0 j Nx j i i j trong đó fi fxt( ij , ). Chi tiết, với j 0, , Nt 1 ta có uj 1 u j uj 2 uu j j 1 1  0 1 2 0 i 1,2 f1 , t x uj 1 u j uj 2 uu j j 2 2  1 2 3 0 i 2,2 f2 , t x (13) jj 1 j jj uNN 1 u 1 uu N 2 2 NN 1 u i N 1,x x  x x x f 0 . x t 2 x 2 Do điều kiện biên, uj h( t ) và uj k( t ), hệ phƣơng trình (13) đƣợc viết lại nhƣ sau 0 j Nx j © 2020 Trƣờng Đại học Công nghiệp Thành phố Hồ Chí Minh
3. TÍNH VỮNG, ỔN ĐỊNH VÀ HỘI TỤ CỦA PHƢƠNG PHÁP SAI PHÂN 107 HỮU HẠN CHO PHƢƠNG TRÌNH NHIỆT ujj 1 u 2 uu jj  1 1  1 2 0 i 1, 2f1 2 h( tj ), t x x uj 1 u j uj 2 uu j j i 2,,2 2  1 2 3 f 0 t 2x 2 (14) jj 1 j j uN 1 u N 1uN 2 2u N 1  i N 1,x x  x x fj k(t ) . x t 2 x2 2 x j Sắp xếp lại các số hạng trong hệ phƣơng trình (14) ta đƣợc hệ tƣơng đƣơng nhƣ sau  t  t jj 1 jj j i1, u1u 12( 2 uu 1 2 ) f 1 t 2 ht (j ), x x jj 1  t jjjj i 2, u u(), uu 2 u f t 2 2 2 x 1 2 3 2 (15) jj 1  t j j j  t iNu xNN1, 1 u 1( u NNN 2 2u 1 ) f 1 t k( t j ). x x 2 x x x x 2 x  t Để đơn giản trong trình bày, ta đặt r . Hệ phƣơng trình (15) đƣợc viết lại nhƣ sau 2 x  t j 1 j jjj i1, u1u 1 2 ruru 1 2 ft 1 2 ht(j ), x j 1 jj jjj i 2, u2 ur 2 ururu 1 2 2 3 ft 2 , i 3,ujj1 ur ururu jjjj 2 ft , 3 3 2 3 4 2 (16) jjj 1 j j j iN x 2, u NN 2 ur 2 u N 3 2 ruru NNN 2 1 f 1 t, x x x j 1 j j j j  t iN xNN 1, u 1 u 1 r u N 2 2ruNN 1 f 1 t kt ( j ). x x x x x 2 x j 1 Vậy, ta có đƣợc ui ,i 1 N 1 đƣợc tính nhƣ sau  t j 1 jjj i1, u1( 2 ruru 1) 1 2ft 1 2 ht (j ), x jj 1 jjj i 2, u2 r u 1 ( 2 ruru 1) 2 3 f 2 t, i 3,uj1 r u j ( 2 ruru 1) jjj ft , 3 2 3 4 2 (17) jj 1 jj j iN xNN2, u 2 r u 3 ( 2 ruru 1) NNN 2 1 f 1 t, x x j 1 j j j  t iNu x1,N 1 r u N 2 ( 2 r 1)uNN 1 f 1 t k ( t j ). x x x x 2 x Đặt các ma trận © 2020 Trƣờng Đại học Công nghiệp Thành phố Hồ Chí Minh
4. 108 TÍNH VỮNG, ỔN ĐỊNH VÀ HỘI TỤ CỦA PHƢƠNG PHÁP SAI PHÂN HỮU HẠN CHO PHƢƠNG TRÌNH NHIỆT j  t ft ht( ) j 1 2 j u x 1 j j f2 t u2 N 1 j j j j x 1 . F , u , trong đó F, u j j f t uN 2 Nx 2 x j j  t u Nx 1 fN 1 t kt( j ) x 2 x và ma trận 2r 1 r 0 0 0 0 r 2 r 1 1 0 0 0 0r 2 r 1 r 0 0 N 1 N 1 A , A x x . (18) 0 0 0 0 0 0 0 0rr 2 1 r 0 0 0 0r 2 r 1 Hệ phƣơng trình (15) tƣơng đƣơng với phƣơng trình ma trận nhƣ sau j 1 jj u Au F ,,j  0,Nt 1. (19) j 1 j 1 j 1 Sau khi giải hệ trên ta đƣợc ui ,i 1 N 1 . Giá trị biên u0 h( tj 1 ) và uN k( t j 1 ) . 3 HỘI TỤ CỦA PHƯƠNG PHÁP Bổ đề 1. Cho ma trận vuông A n n có dạng 3 đƣờng chéo nhƣ sau b c 0   0 a b c 0 0 A       , abc,, . (20) 0 0 a b c 0  0 a b Trị riêng và vector riêng của A tƣơng ứng nhƣ sau 1/2 a 1 j sin c n 1 1/2 a 2  j j sin   jb 2accos ,v j c n 1 , jN 1, , . n 1 1/2 a n j sin c n 1 Chứng minh. Xem [2, 3]. j Định nghĩa 2 (Chuẩn của sai số). Với mọi điểm trên lƣới(xi,t j ) (0,1) (0, T ) , gọi ui là giá trị xấp xỉ của nghiệm chính xác u( xi , t j ) . Sai số Errt()( uuxj (,,,,t ), uux j ( t ), u j u (x t )) j0 0 j 1 1 jNNjx x i) Cố định thời gian t j , chuẩn của sai số theo không gian © 2020 Trƣờng Đại học Công nghiệp Thành phố Hồ Chí Minh
5. TÍNH VỮNG, ỔN ĐỊNH VÀ HỘI TỤ CỦA PHƢƠNG PHÁP SAI PHÂN 109 HỮU HẠN CHO PHƢƠNG TRÌNH NHIỆT Nx 2j 2 ‖Errt() ‖ 2 x ( uu ()x , t ) jL  iji i 1 ii) Chuẩn của sai số theo không gian và thời gian Nt N x 2j 2 Err()2 tx (() u u x , t ) ‖ ‖ L  ii j j 1 i 1 3.1 Tính vững Nhìn chung tại các điểm lƣới (xi , t j ) , nghiệm chính xác không thỏa công thức rời rạc (12), sai số sẽ có dạng nhƣ sau u( x,,,,,t )() u x t u ( x t )2()( ux t u x t ) j ij1 ij  ij 1 ij ij 1 f(x , t ) . (21) i t 2x ij j Định lý 2. Xấp xỉ theo công thức (12) có tính vững, nghĩa là  i sẽ hội tụ về 0 khi x và t tiến đến 0 . Chứng minh. Khai triển chuỗi Taylor của hàm u( x , t ) theo biến t tại điểm t j , ta đƣợc u( x ,t ) u()(xux,, t t )( ux,()t ) t ttt i j () t t 2 , i ijtijj2 j suy ra utt( x i,t j ) 2 u(xij,,,t 1 )()()( u x ijtijjj tu x t t 1 t ) (t jj 1 t ) 2 u( x ,t ) uxux( ,,),t)() t t tt i j ( t 2 ij tij 2 vậy, ta rút ra đƣợc ux()(),,,t uxt u () x t ij 1 ij u (x ,t ) ttij tOt (( )2 ) (22) t t i j 2 Tƣơng tự, khai triển chuỗi Taylor của hàm u( x , t ) theo biến x tại điểm xi , ta đƣợc uxx( x i,t j ) 2 u()()()( x,,,tj u x ijxij t u x t x x i ) ( x x i ) 2 u() x,t u () x , t xxxij(x x )3 xxxxij (x x ) 4 6i2 4 i với x xi 1 , ta đƣợc u( x ,t ) u()()()( x,,,t u x t u x t x x )xx i j ( x x )2 ij 1 ijxijii 12 ii 1 uxxx() x i,t j3 u xxxx () x i, t j 4 ()xii 1 x () x ii 1 x 6 12 (23) u( x ,t ) uxt()(), u x ,t xx i j () 2 i j x i j x2 x u ()x,t u () x , t xxx i j( )3 xxxx i j ( )4 , 6x 24 x và với x xi 1 , ta đƣợc © 2020 Trƣờng Đại học Công nghiệp Thành phố Hồ Chí Minh
6. 110 TÍNH VỮNG, ỔN ĐỊNH VÀ HỘI TỤ CỦA PHƢƠNG PHÁP SAI PHÂN HỮU HẠN CHO PHƢƠNG TRÌNH NHIỆT u( x ,t ) u()()()( x,,,t u x t u x t x x )xx i j ( x x )2 ij 1 ijxijii 12 ii 1 uxxx() x i,t j3 u xxxx () x i, t j 4 ()xii 1 x () x ii 1 x 6 12 (24) u( x ,t ) uxt()(), u x ,t xx i j () 2 i j x i j x2 x u ()x,t u () x , t xxx i j( )3 xxxx i j ( )4 , 6x 24 x Từ (23) và (24) ta đƣợc u()2()() x,,,,ttt u x u x u ( x t ) ij 1 ij ij 1 u() x ,t xxxxij () 2 (25) 2x xx i j12 x Từ (21), (22) và (25) ta suy ra đƣợc u() x,t u () x , t j u() x,,,tttt i j  u () x t xxxx i j () 2 ft () x itij xxij x ij 2 12 u( x ,,t ) u (x t ) tt i j t  xxxx i j ( )2 2 12 x 2 O( t ( x ) ), j trong đó, ta sử dụng giả thiết phƣơng trình ut  u xx fxt , . Ta thấy rằng  i sẽ hội tụ về 0 khi x và t tiến đến 0 . Do đó, xấp xỉ của ta có tính chất vững. 3.2 Tính ổn định  t Để thuận tiện, ta đặt r . Theo bổ đề 1, ma trận (18) có trị riêng và vector riêng tƣơng ứng nhƣ sau 2 x 1 j sin n 1 2  j j sin  2r 2rcos ,v n 1 , jN 1 ,  , . jn 1 j n j sin n 1 Do đó, ta có (Nx 1 ) min N 1 1 2r 2rcos x N x   1 2r 2r cos max 1 Nx Nếu / Nx 1, ta có đƣợc 2 min N 1 1 4r r , x N x 2 max  1 1 r , Nx © 2020 Trƣờng Đại học Công nghiệp Thành phố Hồ Chí Minh
7. TÍNH VỮNG, ỔN ĐỊNH VÀ HỘI TỤ CỦA PHƢƠNG PHÁP SAI PHÂN 111 HỮU HẠN CHO PHƢƠNG TRÌNH NHIỆT Điều kiện hội tụ cho min dẫn đến 2 1 , (26) r 2 2 4 Nx với điều hiện này, (19) sẽ ổn định. 4 VÍ DỤ SỐ Hình 1. Nghiệm xấp xỉ hội tụ về nghiệm chính xác. Phần này, tôi sẽ xét 2 ví dụ cho hai trƣờng hợp: đầu tiên là r thỏa điều kiện (26), và trƣờng hợp r không thỏa điều kiện (26) để ta khảo sát tính ổn định của nghiệm. Ví dụ 1. Ta xét phƣơng trình đạo hàm riêng ut u xx ,(, x t) ( 0,1) ( 0, TT), 0 , với điều kiện ban đầu u( x ,0) sin( x ), x (0,1), và cho biết thêm, bài toán có biên Dirichlet u(0, t ) u (1,), t 0 t (0, T ). 1 Cho trƣớc  , dễ dàng kiểm tra đƣợc bài toán này có nghiệm chính xác 16 © 2020 Trƣờng Đại học Công nghiệp Thành phố Hồ Chí Minh
8. 112 TÍNH VỮNG, ỔN ĐỊNH VÀ HỘI TỤ CỦA PHƢƠNG PHÁP SAI PHÂN HỮU HẠN CHO PHƢƠNG TRÌNH NHIỆT 1 2t u(xt , ) e 4 sin(2 x ) j Mục tiêu bài toán ở đây, ta tìm giá trị xấp xỉ cho u( x , t ) , ui u( x ij , t ) . Trường hợp 1. Chọn t sao cho (26) đƣợc thỏa, nghĩa là 1 t 1 t r suy ra 8 16 2x 2 2 x Ta có kết quả so sánh nghiệm chính xác và nghiệm xấp xỉ tại t 4 và ta chọn t4 2 x . Tại t 4 , tôi so sánh giá trị của nghiệm chính xác u( x ): u ( x ,4) với nghiệm xấp xỉ. Khi Nx càng lớn thì nghiệm xấp xỉ càng gần với nghiệm chính xác xem (Hình 1). Trường hợp 2. Chọn t sao cho điểu kiện (26) bị vi phạm, nghĩa là 1 t 1 t r suy ra 8 16 2x 2 2 x Ta có kết quả so sánh nghiệm chính xác và nghiệm xấp xỉ tại t 4 và ta chọn t16 2 x . Khi điều kiện (26) bị vi phạm, nghiệm xấp xỉ không hội tụ về nghiệm chính xác. Điều này phản ánh đúng kết quả lý thuyết đã trình bày, xem (Hình 2). Hình 2. Nghiệm xấp xỉ không hội tụ về nghiệm chính xác. 5 KẾT LUẬN Trong bài này tôi đã trình bày phƣơng pháp số để giải phƣơng trình nhiệt bằng phƣơng pháp sai phân hữu hạn. Để đảm bảo phƣơng pháp hội tụ, tôi khảo sát tính vững và tính ổn định bằng phân tích phổ của ma trận. Kết quả bài báo chỉ ra rằng phƣơng pháp sai phân hữu hạn trong trƣờng hợp này có tính vững và tính ổn định chỉ đúng với một số cách chọn r . Ví dụ số cũng đƣợc thiết lập để minh họa kết quả lý thuyết. © 2020 Trƣờng Đại học Công nghiệp Thành phố Hồ Chí Minh
9. TÍNH VỮNG, ỔN ĐỊNH VÀ HỘI TỤ CỦA PHƢƠNG PHÁP SAI PHÂN 113 HỮU HẠN CHO PHƢƠNG TRÌNH NHIỆT REFERENCES [1] K. Sayevand, "Convergence and stability analysis of modified backward time centered space approach for non- dimensionalizing parabolic equation," J. Nonlinear Sci., pp. 11-17, 2014. [2] V. FACK, G. VANDEN BERGHE, H.E. DE MEYER , "Some finite difference methods for computing eigenvalues and eigenvectors of special two-point boundary value problems," Journal of Computational and Applied Mathematics , pp. 211-217 , 1987. [3] S. Sajavicius, "On the eigenvalue problems for differential operators with coupled boundary conditions," Nonlinear Analysis: Modelling and Control, vol. 15, p. 493–500, 2010. [4] J.Wang, "A model of competitive stock trading volume," Journal of Political Economy, vol. 102, p. 127–168, 1994. Ngày nhận bài: 11/11/2019 Ngày chấp nhận đăng: 19/03/2020 © 2020 Trƣờng Đại học Công nghiệp Thành phố Hồ Chí Minh