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

17<br /> TẠP CHÍ KHOA HỌC CÔNG NGHỆ GIAO THÔNG VẬN TẢI, SỐ 20 - 08/2016<br /> <br /> <br /> GIẢI SỐ CHO PHƯƠNG TRÌNH ĐẠO HÀM RIÊNG<br /> TỰA TUYẾN TÍNH CẤP 1 HAI BIẾN<br /> NUMERICAL SOLUTION OF QUASI-LINEAR PARTIAL DIFFERENTIAL<br /> EQUATION<br /> Huỳnh Văn Tùng<br /> Khoa Cơ Bản<br /> 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<br /> 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 –<br /> 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<br /> so sánh với nghiệm giải tích thông qua ví dụ mẫu.<br /> 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<br /> bậc 2, phương pháp Runge – Kutta bậc 4, phương pháp lưới.<br /> Abstract: This paper presents the combination of numerical methods, includes: Grid - based,<br /> characteristic curves, quadratic spline interpolation and fourth order Runge – Kutta to find numerical<br /> solutions for quasi - linear PDEs. The results will be compared with analytic solution via examples.<br /> Keywords: PDE, first order, quasi-linear, interpolation, quadratic spline interpolation, fourth<br /> order Runge – Kutta method, grid - based method.<br /> <br /> 1. Giới thiệu cứu các thông tin về nghiệm của chúng.<br /> Phương trình đạo hàm riêng cấp 1 Trong bài báo này, chúng ta áp dụng phối<br /> thường nảy sinh trong vật lý lý thuyết, động hợp các phương pháp số, bao gồm: Phương<br /> lực học (mô tả những chuyển động chính pháp lưới, phương pháp đường đặc trưng,<br /> tắc), cơ học liên tục (để ghi lại sự bảo toàn phương pháp nội suy Spline bậc 2 và phương<br /> khối lượng, mô men lực) và quang học (để pháp Runge – Kutta bậc 4 để giải số cho bài<br /> mô tả sóng),… toán (1) – (3).<br /> Trong bài báo này chúng ta xét phương 2. Phương pháp lưới<br /> trình đạo hàm riêng tựa tuyến tính cấp 1 hai Phương pháp lưới là một trong các<br /> biến có dạng: phương pháp số thông dụng để giải bài toán<br /> u t +c(x,t).u x = F(x,t,u) (1) biên đối với các phương trình đạo hàm riêng.<br /> Ý tưởng của phương pháp lưới được thể hiện<br /> (x,t)  D = (a,b)×(0,T)<br /> như sau: Trong miền biến thiên của các biến<br /> Thỏa mãn điều kiện đầu : độc lập, chúng ta tạo ra một lưới nhờ các<br /> u(x,0) = φ(x), a < x < b (2) đường thẳng song song với hai trục tọa độ.<br /> Và điều kiện biên : Điểm giao nhau của các đường thẳng đó gọi<br /> là các nút lưới (điểm lưới).<br /> u (a,t) = g (t), 0 < t < T<br />  (3) Ứng với mỗi bài toán, chúng ta áp dụng<br /> u (b,t) = h (t), 0 < t < T các phương pháp số thích hợp để tìm giá trị<br /> Trong đó: xấp xỉ của hàm số cần tìm tại các điểm lưới.<br /> u (x,t) : Hàm số cần tìm của hai biến; Trong phần này, chúng ta xét lưới đều<br /> x,t ; t : Biến thời gian; được xác định bởi hai họ đường thẳng song<br /> song với các trục tọa độ:<br /> c,F : Các hàm số cho trước.<br />  x = x i = a+i.h x , i = 0,m<br /> u u  (4)<br /> Các ký hiệu: u t = , ux = .<br /> t x  t = t j = j.h t , j = 0,n<br /> Phương trình đạo hàm riêng phi tuyến Trong đó: h x >0, h t >0 là các số đã cho<br /> nói chung là rất phức tạp, nhưng ta có thể sử<br /> dụng những kỹ thuật khác nhau để nghiên (bước lưới theo trục Οx,Οt ), hình 1.<br /> 18<br /> Journal of Transportation Science and Technology, Vol 20, Aug 2016<br /> <br /> <br /> t xa xb  x / (t) = c  x(t),t <br /> tn  T  x(t ) = x , t  [t j-1 ;t j ] (8)<br />  j i<br /> <br /> Dọc theo một đường đặc trưng γi,j , hàm<br /> tj<br /> Ni , j u(x,t) = u(x(t),t) là hàm của một biến t . Khi<br /> ht<br /> đó ta có:<br /> t2 du u u dx<br /> = + = u t +c(x,t)u x (9)<br /> t1 dt t x dt<br /> hx x Thay (9) vào (1) ta được:<br /> 0 x0  a x1 xi xm  b du<br /> t 0<br /> Lớp thời gian = F  x(t),t,u(x(t),t)  (10)<br /> Hình 1. Sơ đồ lưới.<br /> dt<br /> Giá trị u i,j là nghiệm xấp xỉ của bài toán<br /> Xuất phát từ giá trị của hàm u (x,t) tại<br /> các điểm lưới trên lớp thời gian t = 0 : Cauchy sau tại t = t j :<br /> u i,0 = u(x i ,0) = φ(x i ), i = 0,m  du<br /> (5)  =F  x(t),t,u(x(t),t) <br />  dt (11)<br /> Trên mỗi lớp thời gian t = t j (j = 1,n) kế u(x i,2s , t j-1 )  u iA<br /> tiếp, chúng ta đi tìm các xấp xỉ:<br /> với t [t j-1;t j ]<br /> u (x i ,t j )  u i,j (i = 1,m - 1) (6)<br /> Trong đó: x i,2s là hoành độ của Ai,j-1 .<br /> 3. Phương pháp đường đặc trưng<br /> Lớp thời gian xi ,0  xi<br /> Để khảo sát các đặc tính của hàm u (x,t) tj<br /> u ( x i , t j )  ui , j Ni , j<br /> tại một điểm lưới Ni,j (xi ,t j ) bên trong miền<br /> khảo sát ở lớp thời gian t j (hình 2), người ta tại<br /> N i, j xi ,1<br /> t)<br /> ( x,<br /> tìm một đường cong đặc trưng γi,j nằm trong ht<br /> àm<br /> u xi ,2<br /> rị h  i, j<br /> át A i, j <br /> 1<br /> miền khảo sát nối N i,j với một điểm Ai,j-1 h gi n<br /> địn hâ<br /> c hc<br /> Xá địn<br /> nằm ở lớp thời gian t j-1 đã được khảo sát. xi ,2 s 1 Xá<br /> c<br /> <br /> Ba bài toán quan trọng của phương pháp<br /> Ai , j 1 xi ,2 s Lớp thời gian t j 1 xi<br /> đường đặc trưng cần thực hiện là:<br /> Bài toán 1: Xác định tọa độ chân đường Hình 2. Đường đặc trưng xuất phát từ<br /> Ni , j .<br /> đặc trưng Ai,j-1 . 4. Phương pháp nội suy Spline bậc 2<br /> Bài toán 2: Tìm giá trị nội suy Cho n + 1 mốc nội suy cách đều với<br /> u (xi,2s , t j-1 )  u iA tại chân Ai,j-1 thông qua bước h : (x i ,yi ), i = 0,n .<br /> các giá trị đã biết: ui,j-1 (i = 0,m) . Tìm trên mỗi đoạn [x k-1 ,x k ], k = 1,n một<br /> Bài toán 3: Từ giá trị u iA và đường đặc đa thức bậc hai là Pk (x) sao cho nó đi qua<br /> trưng γi,j , chúng ta tìm lại giá trị u i,j tại nút hai đầu mút và trên hai đoạn kề nhau, hai đa<br /> thức được ghép trơn cấp 1 tại điểm chung<br /> N i,j . (nghĩa là đạo hàm cấp 1 của chúng bằng nhau<br /> Phương trình (1) có họ đặc trưng là: tại điểm đó), hình 3.<br /> dt dx y<br /> Ak ( xk , yk )<br /> = (7)<br /> 1 c(x,t) Pk ( x)<br /> Nếu coi x = x(t) là hàm của biến t thì Pk 1 ( x)<br /> x = x(t) là phương trình của đường đặc trưng x<br /> γi,j và nó là nghiệm của bài toán Cauchy: 0 xk 1 xk xk 1<br /> Hình 3. Hai đa thức được ghép trơn tại Ak .<br />  19<br /> TẠP CHÍ KHOA HỌC CÔNG NGHỆ GIAO THÔNG VẬN TẢI, SỐ 20 - 08/2016<br /> <br /> <br /> Tài liệu [7] đã đề nghị: k1i = h.f(x i ,yi )<br /> <br /> k 2i = h.f  x i +0.5h,yi +0.5k1i <br /> x - x k-1 x - xk<br /> Pk (x) = yk - y k-1<br /> h h (12)  (18)<br /> + a k (x - x k-1 )(x - x k ) k 3i = h.f  x i +0.5h,yi +0.5k 2i <br /> k = h.f  x +h,y +k <br /> Trong đó các hệ số ak được tính theo  4i i i 3i<br /> <br /> công thức: 5.2. Giải bài toán 1 (xác định Ai,j-1)<br /> y - 2yk + yk+1 Áp dụng phương pháp Runge – Kutta bậc<br /> a k+1 + a k = k-1 (13) 4 đối với bài toán Cauchy (8) trên đoạn<br /> h2<br /> (với k = 1,2,...,n - 1 ) [t j-1;t j ] để tìm gần đúng chân đường đặc<br /> Hệ (13) có n - 1 phương trình nên còn trưng Ai,j-1 (xi,2s ;t j-1 ) xuất phát từ điểm nút<br /> thiếu một điều kiện để xác định n hệ số Ni,j (xi ,t j ) . Chia đoạn [t j-1 ,t j ] thành 2s đoạn<br /> a k (k = 1,...,n) . Người ta thường bổ sung<br /> ht<br /> thêm một điều kiện ở đầu dãy dữ liệu: nhỏ bằng nhau với bước h tj = , các điểm<br /> y - y - a .h 2s<br /> P1/ (x 0 ) = a 0 Ûa1 = 1 02 0 (14) chia là: t j,l = t j - l.h tj , l = 0,2s . Các công thức<br /> h<br /> Trong đó a 0 tương ứng là độ dốc của dãy Runge – Kutta bậc 4 xác định x i,2s là:<br /> số liệu ở đầu dãy. Chúng được xấp xỉ bởi đạo  x i,0 = x i<br /> hàm của đa thức nội suy Newton tiến qua <br />  1 (19)<br /> bốn mốc nội suy ở đầu dãy số liệu.  i,l+1 i,l<br /> x = x +  k +2k +2k +k <br /> <br /> 1l 2l 3l 4l<br /> 1 6<br /> 1 1 <br /> a 0 =  Δy0 - Δ 2 y0 + Δ 3 y0  (15) (với l = 0,2s -1 )<br /> h 2 3 <br /> Trong đó: Δyi = yi+1 - yi là sai phân cấp k1l = ht j .c(x i,l ,t j,l )<br /> <br /> 1; Δk yi = Δk-1yi+1 - Δk-1yi là sai phân cấp k . k 2l = ht j .c  x i,l +0.5k1l ,t j,l -0.5h tj <br />  (20)<br /> k 3l = ht j .c  x i,l +0.5k 2l ,t j,l -0.5h tj <br /> Phương pháp nội suy Spline bậc 2 được<br /> áp dụng để giải bài toán 2 trong phương<br /> <br /> pháp đường đặc trưng. k 4l = ht j .c  x i,l +k 3l ,t j,l -h tj <br /> 5. Giải bài toán 1 và bài toán 3<br /> 5.1. Phương pháp Runge - Kutta bậc 4 Khi l = 2s - 1 ta thu được x i,2s , đồng thời<br /> Xét bài toán Cauchy: tại các bước của vòng lặp, ta lưu lại các giá<br />  y / = f(x,y) trị x i,1 ,x i,2 ,..., xi,2s-1 để sử dụng cho bài toán<br />  , x0  x  x (16)<br />  y(x 0 ) = y 0 3.<br /> 5.3. Giải bài toán 3 (tìm ui,j)<br /> Chia đoạn [x 0 ,x] thành n đoạn bằng<br /> Khi chân đường đặc trưng Ai,j-1 được<br /> nhau bởi các điểm chia<br /> xác định, áp dụng phương pháp nội suy<br /> x-x 0<br /> x i = x 0 +i.h ;i = 0,n với bước h = . Cần Spline bậc 2 để tìm xấp xỉ giá trị<br /> n u(xi,2s , t j-1 )  u iA tại Ai,j-1 thông qua các giá<br /> tìm các giá trị xấp xỉ y(x i+1 )  yi+1 ,i = 0,n-1 .<br /> trị đã biết: u 0,j-1 , u1,j-1 ,..., u m-1,j-1 , u m,j-1 .<br /> Xuất phát từ giá trị y0 , hai nhà toán học<br /> Áp dụng phương pháp Runge – Kutta<br /> người Đức là Runge và Kutta đề xuất một<br /> bậc 4 cho bài toán (11) trên đoạn [t j-1;t j ] để<br /> phương pháp tìm các giá trị yi+1 ,i = 0,n-1<br /> tìm gần đúng giá trị u i,j tại Ni,j (xi ,t j ) .<br /> theo công thức sau:<br /> 1 Để sử dụng được các giá trị x i,0 ,x i,1 ,...,<br /> yi+1 = yi +  k1i +2k 2i +2k 3i +k 4i  (17)<br /> 6 x i,2s dọc theo đường cong đặc trưng γi,j , ta<br /> Trong đó: chia đoạn [t j-1;t j ] thành s đoạn nhỏ bằng<br /> 20<br /> Journal of Transportation Science and Technology, Vol 20, Aug 2016<br /> <br /> <br /> ht ht φ(x) = cosx , 0 < x < π<br /> nhau với bước chia: Δt j = =2 = 2h tj <br /> s 2s Đặt: g(t) = (1+ t)e t , 0 < t < 1<br /> bởi các điểm chia: h(t) = (1 + t)e t ,0 < t < 1<br /> z j,r = t j-1 + r.Δt j = t j,2s-2r <br /> (21) Xét các bước lưới tương ứng trục Οx,Οt<br /> = t j - (2s - 2r)h tj , r = 0,s π<br /> là: h x = , h t = 0.1<br /> Ta có: 20<br /> z j,r = t j-1 +r.Δt j = t j,2s-2r Tập các điểm lưới của miền khảo sát là:<br />  (x ,t ):x = i.h x ; t j = j.h t ; <br /> z j,r + 0.5Δt j = t j,2s-2r +h tj = t j,2s-2r-1 G = i j i <br />   i = 0,20 ; j = 0,10 <br /> z j,r +Δt j = t j,2s-2r +2h tj = t j,2s-2r-2<br /> Ký hiệu:<br />  x(z j,r )  x i,2s-2r t<br />  Uexact = (1 + t j ).e j .cos(x i ) là giá trị của<br />   x(t j,2s-2r-1 )  x i,2s-2r-1 (22) nghiệm đúng tại nút lưới Ni,j (xi ,t j ) .<br /> <br />  x(z j,r +Δt j )  x i,2s-2r-2  <br /> Uj = ui,j:i = 0,n là tập các giá trị xấp xỉ<br /> Các công thức Runge – Kutta bậc 4 xác<br /> của hàm u(x,t) tại các nút lưới trên lớp thời<br /> định u i,j :<br /> gian t j .<br /> u i,j0 = u iA<br />  Lập trình số với phần mềm Mathematica,<br />  r+1 r 1 (23) ta có kết quả trong bảng sau trên lớp thời<br /> u i,j = u i,j +  k1r +2k 2r +2k 3r +k 4r <br />  6 gian t 7 .<br /> (với r  0, s  1 ) Bảng 1. So sánh giá trị xấp xỉ và giá trị đúng của<br /> hàm u ( x, t ) trên lớp t7  7h t .<br /> k1r = Δt r F  x i,2s-2r ,z j,r ,u i,jr <br />  Xấp xỉ<br /> k 2r = Δt r F  x i,2s-2r-1 ,z j, r +0.5Δt j ,u i,jr +0.5k1r <br /> STT xi Uexact U7<br />  (24)<br /> k 3r = Δt r F  x i,2s-2r-1 ,z j, r +0.5Δt j ,u i,j +0.5k 2r <br /> r<br /> 0 0.00000 3.42338 3.42338<br /> 1 0.15708 3.38123 3.38270<br /> <br /> k 4r = Δt r F  x i,2s-2r-2 ,z j, r +Δt j ,u i,j +k 3r <br /> r<br /> 2 0.31416 3.25583 3.25667<br /> 3 0.47124 3.05025 3.05072<br /> Khi r = s -1 , ta thu được giá trị cần tìm:<br /> 4 0.62832 2.76957 2.76971<br /> u i,j = u i,js .<br /> 5 0.78540 2.42069 2.42062<br /> 6. Giải số 6 0.94248 2.01221 2.01205<br /> Xét phương trình: 7 1.09956 1.55418 1.55402<br /> u t + c(x,t).u x = F(x,t,u) (25) 8 1.25664 1.05788 1.05761<br /> Điều kiện đầu: 9 1.41372 0.53553 0.53493<br /> u(x,0) = cosx , 0 < x < π (26) 10 1.57077 0.00000 -0.00044<br /> Điều kiện biên: 11 1.72788 -0.53553 -0.53483<br /> 12 1.88496 -1.05788 -1.06084<br /> <br /> u(0, t) = (1+ t) e<br /> t<br /> <br />  ; 0 < t