Giải số cho phương trình đạo hàm riêng tựa tuyến tính cấp 1 hai biến - Huỳnh Văn Tùng

TẠP CHÍ KHOA HỌC CÔNG NGHỆ GIAO THÔNG VẬN TẢI, SỐ 20 - 08/2016 17 GIẢI SỐ CHO PHƯƠNG TRÌNH ĐẠO HÀM RIÊNG TỰA TUYẾN TÍNH CẤP 1 HAI BIẾN NUMERICAL SOLUTION OF QUASI-LINEAR PARTIAL DIFFERENTIAL EQUATION Huỳnh Văn Tùng Khoa Cơ Bản Tóm tắt: Trong bài báo này chúng tôi trình bày phối hợp các phương pháp số, bao gồm: phương pháp lưới, phương pháp đường đặc trưng, phương pháp nội suy Spline bậc 2 và phương pháp Runge – Kutta bậc 4 để giải số cho phương trình đạo hàm riêng tựa tuyến tính cấp 1 hai biến. Kết quả số được so sánh với nghiệm giải tích thông qua ví dụ mẫu. Từ khóa: Phương trình đạo hàm riêng, cấp 1, tựa tuyến tính, phép nội suy, phương pháp Spline bậc 2, phương pháp Runge – Kutta bậc 4, phương pháp lưới. Abstract: This paper presents the combination of numerical methods, includes: Grid - based, characteristic curves, quadratic spline interpolation and fourth order Runge – Kutta to find numerical solutions for quasi - linear PDEs. The results will be compared with analytic solution via examples. Keywords: PDE, first order, quasi-linear, interpolation, quadratic spline interpolation, fourth order Runge – Kutta method, grid - based method. 1. Giới thiệu Phương trình đạo hàm riêng cấp 1 thường nảy sinh trong vật lý lý thuyết, động lực học (mô tả những chuyển động chính tắc), cơ học liên tục (để ghi lại sự bảo toàn khối lượng, mô men lực) và quang học (để mô tả sóng), Trong bài báo này chúng ta xét phương trình đạo hàm riêng tựa tuyến tính cấp 1 hai biến có dạng: t xu +c(x,t).u = F(x,t,u) (1) (x,t) D = (a,b)×(0,T) Thỏa mãn điều kiện đầu : u(x,0) = φ(x), a < x < b (2) Và điều kiện biên : u (a,t) = g (t), 0 < t < T u (b,t) = h (t), 0 < t < T    (3) Trong đó: u (x,t) : Hàm số cần tìm của hai biến; x,t ; t : Biến thời gian; c,F : Các hàm số cho trước. Các ký hiệu: t x u u u = , u = t x     . Phương trình đạo hàm riêng phi tuyến nói chung là rất phức tạp, nhưng ta có thể sử dụng những kỹ thuật khác nhau để nghiên cứu các thông tin về nghiệm của chúng. Trong bài báo này, chúng ta áp dụng phối hợp các phương pháp số, bao gồm: Phương pháp lưới, phương pháp đường đặc trưng, phương pháp nội suy Spline bậc 2 và phương pháp Runge – Kutta bậc 4 để giải số cho bài toán (1) – (3). 2. Phương pháp lưới Phương pháp lưới là một trong các phương pháp số thông dụng để giải bài toán biên đối với các phương trình đạo hàm riêng. Ý tưởng của phương pháp lưới được thể hiện như sau: Trong miền biến thiên của các biến độc lập, chúng ta tạo ra một lưới nhờ các đường thẳng song song với hai trục tọa độ. Điểm giao nhau của các đường thẳng đó gọi là các nút lưới (điểm lưới). Ứng với mỗi bài toán, chúng ta áp dụng các phương pháp số thích hợp để tìm giá trị xấp xỉ của hàm số cần tìm tại các điểm lưới. Trong phần này, chúng ta xét lưới đều được xác định bởi hai họ đường thẳng song song với các trục tọa độ: i x j t x = x = a+i.h , i = 0,m t = t = j.h , j = 0,n    (4) Trong đó: x th >0, h >0 là các số đã cho (bước lưới theo trục Οx,Οt ), hình 1. 18 Journal of Transportation Science and Technology, Vol 20, Aug 2016 xh th 0x a mx b x t 0 1t 2t jt nt T 1x ix x a x b 0t Lớp thời gian ,i jN Hình 1. Sơ đồ lưới. Xuất phát từ giá trị của hàm u (x,t) tại các điểm lưới trên lớp thời gian t = 0 : i,0 i iu = u(x ,0) = φ(x ), i = 0,m (5) Trên mỗi lớp thời gian jt = t (j = 1,n) kế tiếp, chúng ta đi tìm các xấp xỉ: i j i,ju (x ,t ) u (i = 1,m - 1) (6) 3. Phương pháp đường đặc trưng Để khảo sát các đặc tính của hàm u (x,t) tại một điểm lưới i,j i jN (x ,t ) bên trong miền khảo sát ở lớp thời gian jt (hình 2), người ta tìm một đường cong đặc trưng i,jγ nằm trong miền khảo sát nối i,jN với một điểm i,j-1A nằm ở lớp thời gian j-1t đã được khảo sát. Ba bài toán quan trọng của phương pháp đường đặc trưng cần thực hiện là: Bài toán 1: Xác định tọa độ chân đường đặc trưng i,j-1A . Bài toán 2: Tìm giá trị nội suy i,2s j-1u (x , t ) iAu tại chân i,j-1A thông qua các giá trị đã biết: i,j-1u (i = 0,m) . Bài toán 3: Từ giá trị iAu và đường đặc trưng i,jγ , chúng ta tìm lại giá trị i,ju tại nút i,jN . Phương trình (1) có họ đặc trưng là: dt dx = 1 c(x,t) (7) Nếu coi x = x(t) là hàm của biến t thì x = x(t) là phương trình của đường đặc trưng i,jγ và nó là nghiệm của bài toán Cauchy:  / j-1 j j i x (t) = c x(t),t , t [t ;t ] x(t ) = x    (8) Dọc theo một đường đặc trưng i,jγ , hàm u(x,t) = u(x(t),t) là hàm của một biến t . Khi đó ta có: t x du u u dx = + = u +c(x,t)u dt t x dt     (9) Thay (9) vào (1) ta được:   du = F x(t),t,u(x(t),t) dt (10) Giá trị i,ju là nghiệm xấp xỉ của bài toán Cauchy sau tại jt = t :   i,2s j-1 iA du =F x(t),t,u(x(t),t) dt u(x , t ) u      (11) với j-1 jt [t ;t ] Trong đó: i,2sx là hoành độ của i,j-1A . ,2ix ix ,i jN , 1i jA  ,0i ix x ,1ix ,2 1i sx  ,2i sx Xá c đ ịnh ch ân jtLớp thời gian 1jt Lớp thời gian ht ,( , )i j i ju x t u Xá c đ ịnh gi á t rị h àm tạ i ( , ) u x t ,i j ,i jN , 1 i jA  Hình 2. Đường đặc trưng xuất phát từ ,i jN . 4. Phương pháp nội suy Spline bậc 2 Cho n + 1 mốc nội suy cách đều với bước h : i i(x ,y ), i = 0,n . Tìm trên mỗi đoạn k-1 k[x ,x ], k = 1,n một đa thức bậc hai là kP (x) sao cho nó đi qua hai đầu mút và trên hai đoạn kề nhau, hai đa thức được ghép trơn cấp 1 tại điểm chung (nghĩa là đạo hàm cấp 1 của chúng bằng nhau tại điểm đó), hình 3. ( , )k k kA x y kx1kx  1kx  x y ( )kP x 1( )kP x 0 Hình 3. Hai đa thức được ghép trơn tại kA . TẠP CHÍ KHOA HỌC CÔNG NGHỆ GIAO THÔNG VẬN TẢI, SỐ 20 - 08/2016 19 Tài liệu [7] đã đề nghị: k-1 k k k k-1 k k-1 k x - x x - x P (x) = y - y h h + a (x - x )(x - x ) (12) Trong đó các hệ số ka được tính theo công thức: k-1 k k+1 k+1 k 2 y - 2y + y a + a = h (13) (với k = 1,2,...,n - 1) Hệ (13) có n - 1 phương trình nên còn thiếu một điều kiện để xác định n hệ số ka (k = 1,...,n) . Người ta thường bổ sung thêm một điều kiện ở đầu dãy dữ liệu: / 1 0 0 1 0 0 1 2 y - y - a .h P (x ) = a Ûa = h (14) Trong đó 0a tương ứng là độ dốc của dãy số liệu ở đầu dãy. Chúng được xấp xỉ bởi đạo hàm của đa thức nội suy Newton tiến qua bốn mốc nội suy ở đầu dãy số liệu. 2 3 0 0 0 0 1 1 1 a = Δy - Δ y + Δ y h 2 3       (15) Trong đó: i i+1 iΔy = y - y là sai phân cấp 1; k k-1 k-1 i i+1 iΔ y = Δ y - Δ y là sai phân cấp k . Phương pháp nội suy Spline bậc 2 được áp dụng để giải bài toán 2 trong phương pháp đường đặc trưng. 5. Giải bài toán 1 và bài toán 3 5.1. Phương pháp Runge - Kutta bậc 4 Xét bài toán Cauchy: / 0 0 0 y = f(x,y) , x x x y(x ) = y     (16) Chia đoạn 0[x ,x] thành n đoạn bằng nhau bởi các điểm chia i 0x = x +i.h ;i = 0,n với bước 0x-xh = n . Cần tìm các giá trị xấp xỉ i+1 i+1y(x ) y ,i = 0,n-1 . Xuất phát từ giá trị 0y , hai nhà toán học người Đức là Runge và Kutta đề xuất một phương pháp tìm các giá trị i+1y ,i = 0,n-1 theo công thức sau:  i+1 i 1i 2i 3i 4i 1 y = y + k +2k +2k +k 6 (17) Trong đó:       1i i i 2i i i 1i 3i i i 2i 4i i i 3i k = h.f(x ,y ) k = h.f x +0.5h,y +0.5k k = h.f x +0.5h,y +0.5k k = h.f x +h,y +k        (18) 5.2. Giải bài toán 1 (xác định Ai,j-1) Áp dụng phương pháp Runge – Kutta bậc 4 đối với bài toán Cauchy (8) trên đoạn j-1 j[t ;t ] để tìm gần đúng chân đường đặc trưng i,j-1 i,2s j-1A (x ;t ) xuất phát từ điểm nút i,j i jN (x ,t ) . Chia đoạn j-1 j[t ,t ] thành 2s đoạn nhỏ bằng nhau với bước t tj h h = 2s , các điểm chia là: j,l j tjt = t - l.h , l = 0,2s . Các công thức Runge – Kutta bậc 4 xác định i,2sx là:   i,0 i i, l+1 i,l 1l 2l 3l 4l x = x 1 x = x + k +2k +2k +k 6      (19) (với l = 0,2s -1)       1l j i,l j, l 2l j i,l 1l j,l tj 3l j i,l 2l j,l tj 4l j i,l 3l j,l tj k = ht .c(x ,t ) k = ht .c x +0.5k ,t -0.5h k = ht .c x +0.5k ,t -0.5h k = ht .c x +k ,t -h        (20) Khi l = 2s - 1 ta thu được i,2sx , đồng thời tại các bước của vòng lặp, ta lưu lại các giá trị i,1 i,2x ,x ,..., i,2s-1x để sử dụng cho bài toán 3. 5.3. Giải bài toán 3 (tìm ui,j) Khi chân đường đặc trưng i,j-1A được xác định, áp dụng phương pháp nội suy Spline bậc 2 để tìm xấp xỉ giá trị i,2s j-1 iAu(x , t ) u tại i,j-1A thông qua các giá trị đã biết: 0,j-1 1,j-1 m-1,j-1 m,j-1u ,u ,...,u ,u . Áp dụng phương pháp Runge – Kutta bậc 4 cho bài toán (11) trên đoạn j-1 j[t ;t ] để tìm gần đúng giá trị i,ju tại i,j i jN (x ,t ) . Để sử dụng được các giá trị i,0 i,1x ,x ,..., i,2sx dọc theo đường cong đặc trưng i,jγ , ta chia đoạn j-1 j[t ;t ] thành s đoạn nhỏ bằng 20 Journal of Transportation Science and Technology, Vol 20, Aug 2016 nhau với bước chia: t t j tj h h Δt = = 2 = 2h s 2s bởi các điểm chia: j, r j-1 j j,2s-2r j tj z = t + r.Δt = t = t - (2s - 2r)h , r = 0,s (21) Ta có: j, r j-1 j j,2s-2r j, r j j,2s-2r tj j,2s-2r-1 j, r j j,2s-2r tj j,2s-2r-2 z = t +r.Δt = t z + 0.5Δt = t +h = t z +Δt = t +2h = t      j, r i,2s-2r j,2s-2r-1 i,2s-2r-1 j, r j i,2s-2r-2 x(z ) x x(t ) x x(z +Δt ) x        (22) Các công thức Runge – Kutta bậc 4 xác định i,ju :   0 i,j iA r+1 r i,j i,j 1r 2r 3r 4r u = u 1 u = u + k +2k +2k +k 6      (23) (với 0, 1r s  )         r 1r r i,2s-2r j,r i,j r 2r r i,2s-2r-1 j, r j i,j 1r r 3r r i,2s-2r-1 j, r j i,j 2r r 4r r i,2s-2r-2 j, r j i,j 3r k = Δt F x ,z ,u k = Δt F x ,z +0.5Δt ,u +0.5k k = Δt F x ,z +0.5Δt ,u +0.5k k = Δt F x ,z +Δt ,u +k        (24) Khi r = s -1, ta thu được giá trị cần tìm: s i,j i,ju = u . 6. Giải số Xét phương trình: t xu + c(x,t).u = F(x,t,u) (25) Điều kiện đầu: u(x,0) = cosx , 0 < x < π (26) Điều kiện biên: t t u(0, t) = (1+ t)e ; 0 < t <1 u(p, t) = -(1+ t)e    (27) Trong đó: -te + sinx c(x,t) = 1 + t -t 2 t 2 e F(x,t,u) = u + u + (cosx - 1)e - sinx (1+t) Nghiệm chính xác của bài toán (25) – (27) là: tu(x,t) = (1 + t).e .cosx (28) Đặt: t t φ(x) = cosx , 0 < x < π g(t) = (1+ t)e ,0 < t < 1 h(t) = (1 + t)e ,0 < t < 1      Xét các bước lưới tương ứng trục Οx,Οt là: x t π h = , h = 0.1 20 Tập các điểm lưới của miền khảo sát là: i j i x j t(x ,t ):x = i.h ; t = j.h ; G = i = 0,20 ; j = 0,10       Ký hiệu: jt j iUexact = (1 + t ).e .cos(x ) là giá trị của nghiệm đúng tại nút lưới i,j i jN (x ,t ) .  i,jUj = u :i = 0,n là tập các giá trị xấp xỉ của hàm u(x,t) tại các nút lưới trên lớp thời gian jt . Lập trình số với phần mềm Mathematica, ta có kết quả trong bảng sau trên lớp thời gian 7t . Bảng 1. So sánh giá trị xấp xỉ và giá trị đúng của hàm ( , )u x t trên lớp 7 7 t t h . STT ix Uexact Xấp xỉ 7U 0 0.00000 3.42338 3.42338 1 0.15708 3.38123 3.38270 2 0.31416 3.25583 3.25667 3 0.47124 3.05025 3.05072 4 0.62832 2.76957 2.76971 5 0.78540 2.42069 2.42062 6 0.94248 2.01221 2.01205 7 1.09956 1.55418 1.55402 8 1.25664 1.05788 1.05761 9 1.41372 0.53553 0.53493 10 1.57077 0.00000 -0.00044 11 1.72788 -0.53553 -0.53483 12 1.88496 -1.05788 -1.06084 13 2.04204 -1.55418 -1.54922 14 2.19911 -2.01221 -2.02083 15 2.35619 -2.42069 -2.40824 16 2.51327 -2.76957 -2.78701 17 2.67035 -3.05025 -3.02799 18 2.82743 -3.25583 -3.28255 19 2.98451 -3.38123 -3.35183 20 3.14159 -3.42338 -3.42338 Nhận xét: Do xuất phát từ các giá trị đúng trên lớp đầu tiên 0t = t = 0 nên khi j TẠP CHÍ KHOA HỌC CÔNG NGHỆ GIAO THÔNG VẬN TẢI, SỐ 20 - 08/2016 21 càng lớn thì các giá trị xấp xỉ i,ju trên lớp jt = t có thể chứa sai số tương đối lớn. Từ giá trị bảng 1 của i,ju so với giá trị đúng i 7u(x ,t ) trên lớp 7t cho ta thấy rằng phương pháp phối hợp các phương pháp số được trình bày trong bài báo cho kết quả xấp xỉ khá tốt. Hơn nữa, các giá trị i,ju được tính chỉ dựa trên các giá trị i,j-1u {i = 0,m} nên độ phức tạp của thuật toán là tương đối thấp. Hình 4. Giá trị đúng và xấp xỉ trên lớp 7t . Các đồ thị sau so sánh hàm nghiệm đúng tu(x,t) = (1 + t).e .cosx và nghiệm xấp xỉ  U = Uj:j = 0,10 . Hình 5. Đồ thị nghiệm đúng ( , )u x t . Hình 6. Đồ thị mặt lưới xấp xỉ U . 7. Kết luận Việc sử dụng kết hợp các phương pháp số như phần trên đã trình bày sẽ giúp chúng ta giải số nhiều lớp các phương trình đạo hàm riêng với độ chính xác cao. Phương pháp cũng có thể mở rộng cho bài toán số đối với phương trình đạo hàm riêng tựa tuyến tính cấp 1 nhiều biến  Tài liệu tham khảo [1] R.S. Varga, Functional Analysis and Approximation Theory in Numerical Analysis, SIAM, Philadelphia, Pensylvania, 1971. [2] W. Rudin, Principles of Mathematical Analysys, (3rd Ed.), McGrawHill, 1976. [3] L.C. Evans, Partial differential Equations, AMS Press,1998. [4] Phạm Kỳ Anh, Giải tích số, NXB Đại học Quốc gia Hà Nội, 1998. [5] Nguyễn Minh Chương (Chủ biên), Giải tích số, NXB Giáo dục, 2003. [6] Trần Đức Vân, Lý thuyết phương trình vi phân đạo hàm riêng, NXB Đại học Quốc gia Hà Nội, 2005. [7] Huỳnh Văn Tùng, Tính tích phân xác định bằng phương pháp Spline bậc 2, Tạp chí KHCN GTVT, số 3 – 11/2012. Ngày nhận bài: 25/07/2016 Ngày chuyển phản biện: 28/07/2016 Ngày hoàn thành sửa bài: 13/08/2016 Ngày chấp nhận đăng: 20/08/2016