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

 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 
 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 
 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 
 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 
 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 
 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 
 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 
 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 
 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