Phương pháp số giải phương trình truyền nhiệt - Ôn Ngũ Minh

Ôn Ngũ Minh Tạp chí KHOA HỌC & CÔNG NGHỆ 102(02): 155 - 159 155 Ngày nhận bài: 03/11/2012, ngày phản biện: 07/3/2013, ngày duyệt đăng: 26/3/2013 PHƯƠNG PHÁP SỐ GIẢI PHƯƠNG TRÌNH TRUYỀN NHIỆT Ôn Ngũ Minh* Trường Đại học Kỹ thuật Công nghiệp – ĐH Thái Nguyên TÓM TẮT Kỹ thuật số giải phương trình đạo hàm riêng gồm hai phương pháp chủ yếu: sai phân hữu hạn và phần tử hữu hạn. Trong bài viết này, chúng tôi quan tâm tới phương pháp sai phân hữu hạn để giải phương trình dạng parabolic, sử dụng phương trình truyền nhiệt làm ví dụ. Chúng tôi cũng sử dụng ngôn ngữ MATLAB để viết mã cho các giải thuật. Từ khóa: sai phân hữu hạn, phương trình truyền nhiệt, ba đường chéo, phân tích LU. Bài toán truyền nhiệt* Nhiệt độ trong thanh mỏng có thể được mô tả bởi phương trình truyền nhiệt một chiều: ut = auxx, 0 < x < X, 0 < t. Nhiệt độ ban đầu tại mỗi điểm trên thanh được cho bởi: u(x, 0) = f(x), 0 ≤ x ≤ X. Ngoài ra, cần có thêm thông tin tại các điểm mút của thanh. Nếu các mút của thanh được giữ ở nhiệt độ xác định (ví dụ, ngâm mỗi mút trong một môi trường chất lỏng với nhiệt độ dao động theo thời gian), thì các điều kiện biên là: u(0, t) = g1(t), u(X, t) = g2(t), 0 < t. Các trạng thái vật l ý cũng sẽ xác định dạng điều kiện biên khác nhau. Ví dụ, giữ một mút của thanh tương ứng với đạo hàm riêng ux bằng 0 sẽ là ux(0,t) = 0 hoặc ux(X, t) = 0. Khả năng thứ ba là một mút của thanh phải được làm mát (ví dụ, tiếp xúc với không khí lạnh). Điều kiện biên tương ứng liên quan đến sự kết hợp của u và ux tại mút thích hợp. Sự kết hợp của điều kiện ban đầu với các điều kiện biên nói trên là ví dụ điển hình cho phương trình đạo hàm riêng dạng parabolic. * Tel: 0913 351286 Giải phương trình truyền nhiệt Xét phương trình truyền nhiệt một chiều: ut = auxx, 0 ≤ x ≤ X, 0 < t ≤ T, với điều kiện đầu: u(x, 0) = f(x), 0 ≤ x ≤ X, và điều kiện biên: u(0, t) = g1(t), u(X, t) = g2(t), 0 < t ≤ T. Chia đoạn [0, X] thành n đoạn bằng nhau có độ dài h = X/n bởi các điểm chia xi = kh, i = 0, 1, 2, ..., n. Chia đoạn [0, T] thành m phần bằng nhau có độ dài k = T/m bởi các điểm chia tj = jk, j = 0, 1, 2, ..., m. Với mọi i = 0, 1, 2, ..., n và j = 0, 1, 2, ..., m ta ký hiệu: ui,j = u(xi, tj), fi = f(xi), x = X x = 0 Nhiệt độ tại thời điểm t x T Ôn Ngũ Minh Tạp chí KHOA HỌC & CÔNG NGHỆ 102(02): 155 - 159 156 g1j = g1(tj), g2j = g2(tj). Phương pháp thứ nhất (Explicit) Tại bước thứ j theo thời gian, ta xấp xỉ ut bởi công thức sai phân tiến và uxx bởi công thức sai phân trung tâm: ut = ( ), 1 ,1 i j i ju uk + − , uxx = ( )1, , 1,21 2i j i j i ju u uh − +− + . Khi đó ( ) ( ), 1 , 1, , 1,21 2i j i j i j i j i jau u u u uk h+ − +− = − + . Đặt r = 2 ak h ta nhận được: ui,j+1 = rui-1,j + (1 – 2r)ui,j + rui+1,j, Phương pháp này đơn giản vì cho phép tính ui,j+1 thông qua ba giá trị của bước trước theo thời gian. Tuy nhiên, để đảm bảo sự ổn định của nghiệm, ta cần chọn h và k sao cho 0 ≤ r ≤ 0.5. Để cài đặt trên MATLAB, ta tổ chức dữ liệu vào một ma trận u gồm m+1 hàng và n+1 cột. Trong đó các giá trị của g1(t) và g2(t) được lưu tương ứng vào cột thứ nhất và cột cuối cùng (không kể vị trí đầu tiên trong cột). Các giá trị của f(x) được lưu tại hàng đầu tiên. Các giá trị ui,j được lưu ở hàng j+1 cột i+1, với i = 1, 2, ..., n – 1, j = 1, 2, ..., m. Như vậy, mỗi giá trị ui,j sẽ được tính dựa vào ba giá trị của hàng trước mà kề với nó. Hàm heat_explicit() sau đây cho phép giải phương trình truyền nhiệt. Trong đó f, g1, g2 là các biểu thức tương ứng với các hàm f(x), g1(t) và g2(t). Hàm trả lại ma trận nghiệm u, véc tơ hàng x và véc tơ hàng t. function [u,x,t] = heat_explicit(f, g1, g2, X, T, n, m, a) %solve heat equation u_t = au_xx % for 0 <= x <= X, 0 < t <= T %by Explicit method %Initial conditions: u(x,0) = f(x) for 0 <= x <= X %Boundary conditions: u(0,t) = g1(t), u(X,t) = g2(t) %n = number of subintervals for x %m = number of subintervals for t h = X/n; k = T/m; r = a*k/h^2; s = 1 - 2*r; x = [0:h:(n-1)*h,X]; t = [0:k:(m-1)*k,T]; u = zeros(m+1, n+1); for j = 1:m+1 u(j,1)= g1(t(j)); u(j,n+1) = g2(t(j)); end for i = 1:n+1 u(1,i) = f(x(i)); end for j = 2:m+1 u(j,2:n) = r*u(j-1,1:n-1) +s*u(j-1,2:n)+r*u(j-1,3:n+1); end Với bài toán cụ thể: ut = uxx, u(x, 0) = x4, 0 ≤ x ≤ 1, u(0, t) = 0, u(1, t) = 1, 0 < t ≤ 1, lời gọi hàm có dạng như sau (m = 50, n = 5, khi đó h = 0.2, k = 0.02 và r = 0.5): >>[u,x,t]=heat_explicit(@(x)x^4,@(t)0,@(t) 1,1,1,5,50,1); t tj+ xi- x xi+ r = 0.5 cho nghiệm ổn định Ôn Ngũ Minh Tạp chí KHOA HỌC & CÔNG NGHỆ 102(02): 155 - 159 157 f(x0) f(x1) f(x2) ... f(xn-1) f(X) g1(t1) u1,1 u2,1 ... un-1,1 g2(t1) g1(t2) u1,2 u2,2 ... un-1,2 g2(t2) ... ... ... ... ... ... g1(tm-1) u1,m-1 u2,m-1 ... un-1,m-1 g2(tm-1) g1(T) u1,m u2,m ... un-1,m g2(T) Hai hình trên biểu thị nhiệt độ của thanh ứng với r = 0.5 và r = 1.0. Trục chiều dài của thanh hướng từ trái qua phải, trục thời gian hướng từ trước ra sau. Qua đó ta thấy được nhược điểm của phương pháp này. Phương pháp thứ hai (Implicit) Bây giờ ta thay đạo hàm riêng cấp hai theo x bởi công thức sai phân trung tâm tại bước thứ j+1 theo thời gian, uxx = ( )1, , 1,21 2i j i j i ju u uh − +− + . Khi đó ( ) ( ), 1 , 1, 1 , 1 1, 121 2i j i j i j i j i jau u u u uk h+ − + + + +− = − + Đặt r = 2 ak h ta nhận được: –rui-1,j+1 + (1 + 2r)ui,j+1 – rui+1,j+1 = ui,j, i = 1, 2, ..., n – 1. Trong phương pháp này, để tính ui, j+1 ta chỉ sử dụng được một giá trị của bước trước và phải sử dụng hai giá trị chưa biết trong cùng bước. Do đó chúng ta phải giải một hệ phương trình đại số. Để rõ hơn, ta xét bài toán ut = uxx, u(x, 0) = x4, 0 ≤ x ≤ 1, u(0, t) = 0, u(1, t) = 1, 0 < t ≤ 0.2 trong trường hợp m = 5 và n = 5. Khi đó h = 0.2, k = 0.04 nên r = 1, phương trình sai phân được viết lại dưới dạng: –ui-1,j+1 + 3ui,j+1 – ui+1,j+1 = ui,j, i = 1, 2, ..., n – 1. Ta có [x1, x2, x3, x4] = [0.2, 0.4, 0.6, 0,8] nên u1,0 = 0.0016, u2,0 = 0.0256, u3,0 = 0.1296, u4,0 = 0.4096. Vì u(0, t) = 0 và u(1, t) = 1 nên u0,1 = u0,2 = ... = u0,5 = 0, u5,1 = u5,2 = ... = u5,5 = 1. Tại bước j = 0, lần lượt lấy i = 1, 2, 3, 4 ta nhận được hệ: 3u1,1 – u2,1 = u1,0 + u0,1 = 0.0016 + 0.0, –u1,1 + 3u2,1 – u3,1 = u2,0 = 0.0256 – u2,1 + 3u3,1 – u4,1 = u3,0 = 0.1296 – u3,1 + 3u4,1 = u4,0 + u5,1 = = 0.4096 + 1.0 Nhận thấy đây là hệ dạng ba đường chéo (Tridiagonal), nên có thể sử dụng phương pháp phân tích LU để giải hệ này. Để tiện theo dỗi, ta xét ma trận T dạng ba đường chéo cấp bốn với phân tích T = LU: t tj+1 xi-1 xi xi+1 r = 1.0 cho nghiệm không ổn định Ôn Ngũ Minh Tạp chí KHOA HỌC & CÔNG NGHỆ 102(02): 155 - 159 158 L = 2 3 4 1 0 0 0 1 0 0 0 1 0 0 0 1 Lb Lb Lb             , U = 1 1 2 2 3 3 4 0 0 0 0 0 0 0 0 0 Ud U Ud U Ud U Ud             , T = 1 1 2 2 2 3 3 3 4 4 0 0 0 0 0 0 d a b d a b d a b d             Dễ thấy rằng các phần tử của các ma trận L và U thỏa mãn: Ui = ai, Ud1 = d1, Lb2 = b1/Ud1, Ud2 = d2 – Lb2a1, Lb3 = b3/Ud2, Ud3 = d3 – Lb3a2, Lb4 = b4/Ud3, Ud4 = d4 – Lb4a3, Tổng quát lên với T là ma trận vuông cấp n: Ui = ai, i = 1, 2, ..., n-1 và Ud1 = d1. Sau đó, với mỗi i từ 2 đến n ta tính: Lbi = bi/Udi-1, Udi = di – Lbiai-1. Hàm lu_tridiag() sau đây trả lại hai véc tơ Ud và Lb: function [Ud Lb] = lu_tridiag(a, d, b) %LU factorization of a tridiagonal matrix T, T = LU %Input: % a vector of elements above main diagonal, a(end) = 0 % d diagonal of matrix T % b vector of elements below main diagonal, b(1) = 0 %Output: % Ud diagonal of upper bidiagonal matrix U, that upper diagonal is a % Lb Lower diagonal of lower bidiagonal matrix L, % that 1s on main diagonal, Lb(1) = 0 n = length(d); Ud(1) = d(1); for i = 2:n Lb(i) = b(i)/Ud(i-1); Ud(i) = d(i) - Lb(i)*a(i-1); end Như vậy, thay vì giải trực tiếp Tx = r, ta sẽ giải hệ tương đương LUx = r. Hệ này tương đương với hai hệ Ly = r và Ux = y. Hàm lu_triadiag_solve() sau đây thực hiện điều đó. Ở đây a là véc tơ nằm ngay trên đường chéo chính của U, d là véc tơ đường chéo chính của U, còn b là véc tơ dưới đường chéo chính của L. function x = lu_tridiag_solve(a, d, b, r) %Function to solve the equation Tx = r, %where T is tridiagonal matrix %LU factorization of T is expressed as % Lower bidiagonal matrix: diagonal is 1s, % lower diagonal is b, b(1) = 0 % Upper bidiagonal matrix: diagonal is d, % upper diagonal is a, a(n) = 0 %Right hand side vector is r n = length(d); y(1) = r(1); for i = 2:n y(i) = r(i) - b(i)*y(i-1); end x(n) = y(n)/d(n); for i = n-1:-1:1 x(i) = (y(i)-a(i)*x(i+1))/d(i); end Để cài đặt trên MATLAB, ta tổ chức dữ liệu như trong phương pháp trước. Nhận thấy rằng ma trận của hệ tại mỗi bước có dạng ba đường chéo và không đổi, chỉ thay đổi vế phải, vì vậy việc phân tích LU là một lần duy nhất. Hàm heat_implicit() sau đây có các tham số vào và ra giống như hàm heat_explicit(). function [u x t] = heat_implicit(f,g1,g2,X,T,n,m,a) %solve heat equation u_t = au_xx % for 0 <= x <= X, 0 < t <= T %by Implicit method Ôn Ngũ Minh Tạp chí KHOA HỌC & CÔNG NGHỆ 102(02): 155 - 159 159 %Initial conditions: u(x,0) = f(x) for 0 <= x <= X %Boundary conditions: u(0,t) = g1(t), u(X,t) = g2(t) %n = number of subintervals for x %m = number of subintervals for t h = X/n;k = T/m; r = a*k/h^2; x = [0:h:h*(n-1),X]; t = [0:k:k*(m-1),T]; u = zeros(m+1,n+1); for j = 1:m+1 u(j,1) = g1(t(j)); u(j,n+1) = g2(t(j)); end for i = 1:n+1 u(1,i) = f(x(i)); end d = (1+2*r)*ones(1,n-1); aa = -r*ones(1,n-1); aa(n-1) = 0; b = -r*ones(1,n-1); b(1) = 0; [Ud, Lb] = lu_tridiag(aa,d,b); for j = 2:m+1 rh(1) = u(j-1,2) + r*u(j,1); rh(2:n-2) = u(j-1,3:n-1); rh(n-1) = u(j-1,n) + r*u(j,n+1) ; u(j,2:n) = lu_tridiag_solve(aa,Ud,Lb,rh); end Với bài toán: ut = uxx, u(x, 0) = x4, 0 ≤ x ≤ 1, u(0, t) = 0, u(1, t) = 1, 0 < t ≤ 1, lời gọi hàm có dạng như sau (m = 50, n = 50: >> [u,x,t] = heat_implicit(... @(x) x^4, @(t) 0, @(t) 1, 1, 1, 50, 50, 1); Lời gọi sau ứng với m = n = 10: >> [u,x,t] = heat_implicit(... @(x) x^4, @(t) 0, @(t) 1, 1, 1, 10, 10, 1); Kết luận Ta thấy rằng phương pháp thứ hai luôn cho nghiệm ổn định, bất kể ta chọn lưới như thế nào. Tuy nhiên khối lượng tính toán là không nhỏ. Trong tương lai, chúng tôi tiếp tục sử dụng ngôn ngữ MATLAB để viết mã cho các giải thuật sử dụng phương pháp sai phân hữu hạn để giải các bài toán về phương trình đạo hàm riêng như bài toán elliptic và bài toán Poisson. TÀI LIỆU THAM KHẢO [1]. Laurene V. Fausett, Applied Numerical Analysis Using MATLAB, Prentic Hall International, Inc. 1999. SUMMARY NUMERICAL METHODS FOR SOLVING HEAT TRANSFER EQUATIONS On Ngu Minh* College of Technology - TNU * Tel: 0913 351286 Ôn Ngũ Minh Tạp chí KHOA HỌC & CÔNG NGHỆ 102(02): 155 - 159 160 Numerical techniques for solving partial differential equations are primarily of two types: finite- difference method and finite-element method. In this paper, we investigate finite-difference methods for parabolic equations, using the heat equation as an example. We also use the MATLAB language to write code for the algorithms. Keywords: Finite-difference method, Heat equation, Tridiagonal, LU factorization. Ngày nhận bài:04/1/2013, ngày phản biện:26/2/2013, ngày duyệt đăng:26/3/2013