CHƯƠNG 9: PHƯƠNG TRÌNH VI PHÂN ĐẠO HÀM RIÊNG

Chia sẻ: Hikaru Hikaru | Ngày: | Loại File: PDF | Số trang:35

0
537
lượt xem
164
download

CHƯƠNG 9: PHƯƠNG TRÌNH VI PHÂN ĐẠO HÀM RIÊNG

Mô tả tài liệu
  Download Vui lòng tải xuống để xem tài liệu đầy đủ

TÀI LIỆU THAM KHẢO VỀ PHƯƠNG TRÌNH VI PHÂN ĐẠO HÀM RIÊNG

Chủ đề:
Lưu

Nội dung Text: CHƯƠNG 9: PHƯƠNG TRÌNH VI PHÂN ĐẠO HÀM RIÊNG

  1. CHƯƠNG 9: PHƯƠNG TRÌNH VI PHÂN ĐẠO HÀM RIÊNG §1. KHÁI NIỆM CHUNG    Phương trình vi phân đạo hàm riêng(PDE) là một lớp các phương trình  vi  phân  có  số  biến  độc  lập  lớn  hơn  1.  Trong  chương  này  ta  sẽ  khảo  sát  các  phương  trình  vi  phân  đạo  hàm  riêng  cấp  2  với  hai  biến  độc  lập  x  và  y,  có  dạng tổng quát:  ∂2u ∂ 2u ∂ 2u ⎛ ∂u ∂u ⎞   A(x,y) 2 + B(x,y) + C(x, y) 2 = f ⎜ x,y,u, , ⎟     (1)  ∂x ∂x∂y ∂y ⎝ ∂x ∂y ⎠ với xo ≤ x ≤ xf, yo ≤ y ≤ yf và các điều kiện biên:    u(x,y o ) = b yo (x)     u(x,y f ) = b yf (x)    u(xo ,y) = bxo (y) u(xf ,y) = bxf (y) (2)  Các PDE được phân thành 3 loại:    PDE elliptic:  B2 − 4AC < 0     PDE parabolic:  B2 − 4AC = 0     PDE hyperbolic:  B2 − 4AC > 0   Các phương trình này gắn một cách tương ứng với trạng thái cân bằng, trạng  thái truyền nhiệt, hệ thống dao động    §2. PHƯƠNG TRÌNH ELLIPTIC    Ta xét phương trình Helmholz:  ∂ 2 u(x,y) ∂ 2 u(x,y)   ∇ 2 u(x,y) + g(x, y) = + + g(x,y)u(x,y) = f(x,y)   (1)  ∂x 2 ∂y 2 trên miền  D = {(x,y) : xo ≤ x ≤ xf ,y o ≤ y ≤ y f }  với điều kiện biên dạng:  u(x,y o ) = b yo (x) u(x, y f ) = b yf (x)           (2)  u(xo ,y) = b xo (y) u(xf ,y) = b xf (y) Phương  trình  (1)  được  gọi  là  phương  trình  Poisson  nếu  g(x,  y)  =  0  và  gọi  là  phương trình Laplace nếu g(x, y) = 0 và f(x, y) = 0. Để dùng phương pháp sai  phân ta chia miền thành Mx đoạn, mỗi đoạn dài ∆x = (xf ‐ xo)/Mx dọc theo trục  x và thành My đoạn, mỗi đoạn dài ∆y = (yf ‐ yo)/My dọc theo trục y và thay đạo  hàm bậc 2 bằng xấp xỉ 3 điểm:  ∂ 2 u(x,y) u i ,j+1 − 2u i ,j + u i ,j−1   ≅  với xj = xo + j∆x, yj = yo + j∆y  (3.a)  ∂x 2 x j ,yi ∆x 2 403
  2. ∂ 2 u(x,y) u i+1,j − 2u i ,j + u i−1,j ≅  với ui,j = u(xj, yi)      (3.b)  ∂y 2 x ,y ∆x 2 j i Như vậy tại mỗi điểm bên trong (xj, xi) với 1 ≤ i ≤ My ‐ 1 và 1 ≤ j ≤ Mx ‐ ậit nhận  được phương trình sai phân:  u i ,j+1 − 2u i ,j + u i ,j−1 u i+1,j − 2u i ,j + u i−1,j   + = g i ,ju i ,j = fi ,j       (4)  ∆x 2 ∆y 2 Trong đó:    ui,j = u(xj, yi)  fi,j = f(xj, yi)    gi,j = g(xj, yi)  Các phương trình này sắp xếp lại theo cách nào đó thành hệ phương trình với  { (My ‐ 1)(Mx ‐ 1) biến  u1,1 ,u1,2 ,...,u1,Mx −1 ,u 2 ,1 ,...,u 2 ,Mx −1 ,...,u My −1,2 ,...,u My −1,Mx −1 . Để  } dễ dàng ta viết lại phương trình và điều kiện biên dưới dạng:    u i ,j = ry (u i ,j+1 + u i ,j−1 ) + rx (u i+1,j + u i−1,j ) + rxy (g i ,ju i ,j − fi ,j )       (5a)  u i ,o = b xo (y i ) u i ,Mx = b xf (y i ) u o ,j = b yo (x j ) u My ,j = b yf (x j )     (5b)  Trong đó:  ∆y 2 ∆x 2 ∆x 2 ∆y 2   ry =        rx =   rxy =   (6)  2( ∆x 2 + ∆y 2 ) 2( ∆x 2 + ∆y 2 ) 2( ∆x 2 + ∆y 2 ) Bây giờ ta khảo sát tiếp các dạng điều kiên biên. Các bài toán PDE có 2 loại  điều kiện biên: điều kiên biên Neumann và điều kiên biên Dirichlet. Điều kiện  biên Neumann mô tả bằng:  ∂u(x,y)   = b′xo (y)                 (7)  ∂x x=xo Thay đạo hàm bậc 1 ở biên trái (x = xo) bằng xấp xỉ 3 điểm:  u i ,1 − u i ,−1   = b′xo (y i )   u i ,−1 ≈ u i ,1 − 2b′xo (y i )∆x   i = 1, 2,..., My‐1  (8)  2 ∆x Thay thế ràng buộc này vào (5a) ở các điểm biên ta có:    u i ,0 = ry (u i ,1 + u i ,−1 ) + rx (u i+1,0 + u i−1,0 ) + rxy (g i ,0 u i ,0 − fi ,0 )             = ry ⎡ u i ,1 + u i ,1 − 2b′x0 (y i )∆x ⎤ + rx (u i+1,0 + u i−1,0 ) + rxy (g i ,0 u i ,0 − fi ,0 )   ⎣ ⎦          = 2ry u i ,1 + rx (u i+1,0 + u i−1,0 ) + rxy ⎡g i ,0 u i ,0 − fi ,0 − 2b′x0 (y i )∆x ⎤    ⎣ ⎦ (9)  Nếu điều kiên biên trên biên dưới (y = yo) cũng là kiểu Neumann ta sẽ viết các  phương trình tương tự với j = 1, 2,...,Mx‐1:    u 0 ,j = 2rx u1,j + ry (u 0 ,j+1 + u 0 ,j−1 ) + rxy ⎡g 0 ,j u 0 ,j − f0 ,j − 2b′y0 (x j )∆y ⎤    ⎣ ⎦ (10)  và bổ sung cho góc dưới trái(xo, yo):  404
  3. ⎡ b′x (y 0 ) b′y (x0 ) ⎤   u 0 ,0 = 2(ry u 0 ,1 + rx u1,0 ) + rxy ⎢g 0 ,0 u 0 ,0 − f0 ,0 − 2 0 +2 0 ⎥  (11)  ⎣ ∆x ∆y ⎦ Điều kiện biên Dirichlet cho giá trị hàm trên biên nên có thể thay trực tiếp vào  phương trình. Ta có thể lấy giá trị trung bình của các giá trị biên làm giá trị  đầu của ui,j. Ta xây dựng hàm poisson() để thực hiện thuật toán này:    function [u, x, y] = poisson(f, g, bx0, bxf, by0, byf, D, Mx, My, tol, maxiter)  % giai a(u_xx + u_yy + g(x,y)u = f(x,y)  % tren mien D = [x0, xf, y0, yf] = {(x,y) |x0 
  4.     for j = 1:Mx          F(i, j) = f(x(j), y(i));           G(i, j) = g(x(j), y(i));      end  end  dx2 = dx*dx;   dy2 = dy*dy;   dxy2 = 2*(dx2 + dy2);  rx = dx2/dxy2;   ry = dy2/dxy2;   rxy = rx*dy2;  for itr = 1:maxiter      for j = 2:Mx          for i = 2:My              u(i, j) = ry*(u(i, j + 1)+u(i,j ‐ 1)) + rx*(u(i + 1,j)+u(i ‐ 1,j))...                      + rxy*(G(i,j)*u(i,j)‐ F(i,j)); %Pt.(5a)          end      end      if itr > 1 & max(max(abs(u ‐ u0))) 
  5. f = inline(ʹ0ʹ,ʹxʹ,ʹyʹ);   g = inline(ʹ0ʹ,ʹxʹ,ʹyʹ);  x0 = 0;   xf = 4;   Mx = 20;   y0 = 0;   yf = 4;   My = 20;  bx0 = inline(ʹexp(y) ‐ cos(y)ʹ,ʹyʹ); %(vd.2a)  bxf = inline(ʹexp(y)*cos(4) ‐ exp(4)*cos(y)ʹ,ʹyʹ); %(vd.2b)  by0 = inline(ʹcos(x) ‐ exp(x)ʹ,ʹxʹ); %(vd.3a)  byf = inline(ʹexp(4)*cos(x) ‐ exp(x)*cos(4)ʹ,ʹxʹ); %(vd.3b)  D = [x0 xf y0 yf];   maxiter = 500;   tol = 1e‐4;  [U, x, y] = poisson(f, g, bx0, bxf, by0, byf, D, Mx, My, tol, maxiter);  clf  mesh(x, y, U)  axis([0 4 0 4 ‐100 100])     §3. PHƯƠNG TRÌNH PARABOLIC  1.  Dạng  phương  trình:  Một  phương  trình  vi  phân  đạo  hàm  riêng  dạng  parabolic là phương trình mô tả sự phân bố nhiệt độ ở điểm x tại thời điểm t  của một thanh:  ∂ 2 u(x,t) ∂u(x,t)   A =                 (1)  ∂x 2 ∂t Để  phương  trình  có  thể  giải  được  ta  phải  cho  điều  kiện  biên  u(0,  t)  =  b0(t),  u(xf ,t) = bxf (t)  và điều kiện đầu u(x, 0) = i0(x)   2. Phương pháp Euler tiến tường minh: Để áp dụng phương pháp sai phân  hữu  hạn,  ta  chia  miên  không  gian  [0,  xf]  thành  M  đoạn,  mỗi  đoạn  dài  ∆x = xf / M  và chia thời gian T thành N phần, mỗi phần là ∆t = T/N. Sau đó ta  thay đạo hàm bậc 2 ở vế trái và đạo hàm bậc ở vế phải của (1) bằng các xấp xỉ  3 điểm và nhạn được:  u k − 2u ik + u ik−1 u ik+1 − u ik   A i +1 =              (2)  ∆x 2 ∆t 407
  6. Công thức này có thể gói gọn vào thuật toán sau, gọi là thuật toán Eulẻ tiến  tường minh:  ∆t   u ik+1 = r(u ik+1 + u ik−1 ) + (1 − 2r)u ik r = A 2         (3)  ∆x   i = 1, 2,...,M‐1  Để tìm điều kiện ổn định của thuât toán, ta thay nghiệm thử:  jiπ   u =λ e     k i k P                 (4)  với P là số nguyên khác zero vào phương trình (3) và có:  jπ jπ λ = r⎜ ⎛ e P + e − P ⎞ + (1 − 2r) = 1 − 2r ⎡1 − cos π ⎤     ⎟ ⎢       (5)  ⎝ ⎠ ⎣ P⎥⎦ Do ta phải có |λ|≤ 1 với bài toán không có nguồn nên điều kiện ổn định là:  ∆t 1   r=A 2 ≤                   (6)  ∆x 2 Ta xây dựng hàm fwdeuler() để thực hiện thuật toán trên    function [u, x, t] = fwdeuler(a, xf, T, it0, bx0, bxf, M, N)  %giai au_xx = u_t voi  0 
  7.     end  end  3.  Phương  pháp  Euler  lùi  ẩn:  Ta  khảo  sát  một  thuật  toán  khác  gọi  là  thuật  toán Euler lùi, ẩn sinh ra do thay thế lùi xấp xỉ đạo hàm đối với đạo hàm bậc  1 trên vế phải của (1):  u ik+1 − 2u ik + u ik−1 u ik − u ik −1 A =              (7)  ∆x 2 ∆t ∆t   −ru ik−1 + (1 + 2r)u ik − ru ik+1 = u ik −1 r = A 2         (8)  ∆x k k Nếu  các  giá  trị  u 0   và  u M   ở  cả  hai  đầu  đã  cho  trước  từ  điều  kiện  biên  kiểu  Dirichlet nên phương trình (8) đưa tới hệ phương trình:  k −1 ⎡1 + 2r 0 ⎤ ⎡ u1 ⎤ ⎡ u1 + ru 0 ⎤ k k −r 0 0 0 ⎢ ⎥⎢ k ⎥ ⎢ k −1 ⎥ ⎢ −r 1 + 2r −r 0 0 0 ⎥ ⎢ u2 ⎥ ⎢ u2 ⎥ ⎢ 0 −r 1 + 2r −r 0 0 ⎥ ⎢ u k ⎥ ⎢ u k−1 ⎥   ⎢ ⎥ ⎢ 3 ⎥=⎢ 3 ⎥ (9)  ⎢ M M M M M M ⎥⎢ M ⎥ ⎢ M ⎥ ⎢ 0 0 0 L 1 + 2r ⎥⎢ k ⎥ ⎢ −r ⎥ ⎢ u ⎥ ⎢ u k−1 ⎥ ⎢ M−2 M −2 ⎥ ⎢ 0 0 0 L −r 1 + 2r ⎥ ⎢ u k ⎥ ⎢ u k−1 + ru k ⎥ ⎣ ⎦ ⎣ M−1 ⎦ ⎣ M −1 M⎦ ∂u Điều kiện biên Neumann  = b′ (t)  được đưa vào phương trình bằng cách  ∂x x=0 0 xấp xỉ:  u1 − u k1 k   − = b′ (k)                   (10)  2 ∆x 0 k và ghép nó với phương trình có ẩn u 0 :    −ru k1 + (1 + 2r)u 0 − ru1 = u 0 −1     − k k k           (11)  để có được phương trình:    (1 + 2r)u 0 − 2ru1 = u 0 −1 − 2rb′ (k)∆x     k k k 0         (12)  Kết quả ta có được hệ phương trình:  409
  8.   ⎡1 + 2r −r 0 0 L 0 0 ⎤ ⎡ u 0 ⎤ ⎡ u 0 − 2rb′ (k)∆x ⎤ k k 0 ⎢ ⎥⎢ k⎥ ⎢ ⎥ ⎢ −r 1 + 2r −r k −1 0 L 0 0 ⎥ ⎢ u1 ⎥ ⎢ u1 ⎥ ⎢ 0 −r 1 + 2r −r L 0 0 ⎥ ⎢u2 ⎥ ⎢ k u2k −1 ⎥ ⎢ ⎥⎢ k⎥ ⎢ ⎥ ⎢ 0 0 −r 1 + 2r L 0 0 ⎥ ⎢u3 ⎥ = ⎢ u k −1 3 ⎥ (13  ⎢ M M M M M −r 0 ⎥⎢ M ⎥ ⎢ M ⎥ ⎢ ⎥⎢ k⎥ ⎢ ⎥ ⎢ 0 0 0 0 L 1 + 2r −r ⎥ ⎢ u 0 ⎥ ⎢ u M−−2 k 1 ⎥ ⎢ 0 0 0 0 L −r 1 + 2r ⎥ ⎢ u 0 ⎦ ⎢ u M −1 + ru M ⎥ k⎥ k −1 k ⎣ ⎦⎣ ⎣ ⎦ Điểu kiện ổn định của nghiệm là:  jπ jπ − 1   −re P + (1 + 2r) − re = P   λ 1 hay:  λ =   λ ≤ 1             (14)  1 + 2r ⎢⎡1 − cos π ⎤ ⎣ P⎥ ⎦ Ta xây dựng hàm backeuler() để thực hiện thuật toán này:    function [u, x, t] = backeuler(a, xf, T, it0, bx0, bxf, M, N)  %Giai au_xx = u_t voi 0 
  9.     A(i, i) = r2; %Pt.(9)      if i > 1          A(i ‐ 1, i) = ‐r;           A(i, i ‐ 1) = ‐r; end  end  for k = 2:N + 1      b = [r*u(1, k); zeros(M ‐ 3, 1); r*u(M + 1, k)] + u(2:M, k ‐ 1); %Pt.(9)      u(2:M, k) = trid(A, b);  end    4. Phương pháp  Crank ‐ Nicholson: Trong (7), xấp xỉ đạo hàm ở vế trái lấy ở  thời điểm k, trong khi xấp xỉ đạo hàm ở vế phải. Để cải thiện, ta lấy đạo hàm  ở vế trái là trong bình của xấp xỉ đạo hàm tại hai điểm là k và k+1 và có:  A ⎛ u ik++11 − 2u ik+1 + u ik−+11 u ik+1 − 2u ik + u ik−1 ⎞ u ik+1 − u ik + ⎟=       (15)  2⎜ ⎝ ∆x 2 ∆x 2 ⎠ ∆t và nhận được phương pháp Crank ‐ Nicholson:  ∆t   −ru ik++11 + (1 + 2r)u ik+1 − ru ik−+11 = ru ik+1 + (1 − 2r)u ik + ru ik−1 r = A 2   (16)  ∆x Với điều kiện biên Dirichlet tại x0 và điều kiện biên Neumann tại xM ta có hệ  phương trình:  k +1 ⎡ 2(1 + r) −r 0 0 0 0 ⎤ ⎡ u1 ⎤ ⎢ ⎥ ⎢ k +1 ⎥ ⎢ −r 2(1 + r) −r 0 0 0 ⎥ ⎢ u2 ⎥ ⎢ 0 −r 2(1 + r) −r 0 0 ⎥ ⎢ u k +1 ⎥    ⎢ ⎥⎢ 3 ⎥  ⎢ M M M M M M ⎥⎢ M ⎥ ⎢ 0 0 0 L 2(1 + r) ⎥⎢ − r ⎥ ⎢ u k +1 ⎥ ⎥ ⎢ M −1 ⎢ 0 0 0 L −r 2(1 + r) ⎥ ⎢ u k +1 ⎥ ⎣ ⎦⎣ M ⎦ ⎡ 2(1 − r) 0 ⎤ ⎡ u1 ⎤ k r 0 0 0 ⎢ ⎥⎢ k ⎥ ⎢ r 2(1 − r) r 0 0 0 ⎥ ⎢ u2 ⎥ ⎢ 0 r 2(1 − r) r 0 0 ⎥ ⎢ uk ⎥          = ⎢ ⎥⎢ 3 ⎥  ⎢ M M M M M M ⎥⎢ M ⎥ ⎢ 0 ⎢ ⎥ ⎢ 0 0 L 2(1 − r) r ⎥ ⎢uk ⎥ ⎥ M −1 ⎢ 0 0 0 L r 2(1 − r) ⎥ ⎢ u k ⎥ ⎣ ⎦⎣ M ⎦ 411
  10. ⎡ r(u 0 +1 + u 0 ) k k ⎤ ⎢ ⎥ ⎢ 0 ⎥ ⎢ 0 ⎥          + ⎢ ⎥              (17)  ⎢ M ⎥ ⎢ 0 ⎥ ⎢ ⎥ ⎢ 2r [ b′ (k + 1) + b′ (k)]⎥ ⎣ M M ⎦ Điều kiện ổn định được xác định bằng:  ⎡ π ⎤ ⎡ π ⎤   2λ ⎢1 + r ⎛ 1 − cos ⎞ ⎥ = 2 ⎢1 − r ⎛ 1 − cos ⎞ ⎥   ⎜ ⎟ ⎜ ⎟ ⎣ ⎝ P ⎠⎦ ⎣ ⎝ P ⎠⎦ π 1 − r ⎡1 − cos ⎤ ⎢ hay:  λ = ⎣ P⎥ ⎦ λ ≤ 1              (18)  1+ r⎢ ⎡1 − cos π ⎤ ⎣ P⎥ ⎦ Ta xây dựng hàm cranknicholson() để thực hiện thuật toán trên:    function [u, x, t] = cranknicholson(a, xf, T, it0, bx0, bxf, M, N)  %Giai au_xx = u_t voi 0 
  11.     A(i, i) = r2; %Pt.(17)      if i > 1          A(i ‐ 1, i) = ‐r;           A(i, i ‐ 1) = ‐r;       end  end  for k = 2:N + 1      b = [r*u(1, k); zeros(M ‐ 3, 1); r*u(M + 1, k)] ...      + r*(u(1:M ‐ 1, k ‐ 1) + u(3:M + 1, k ‐ 1)) + r1*u(2:M, k ‐ 1);      u(2:M, k) = trid(A,b); %Pt.(17)  end    Để giải phương trình:  ∂ 2 u(x,t) ∂u(x,t)   = 0 ≤ x ≤ 1, 0 ≤ t ≤ 0.1           (vd1)  ∂x 2 ∂t với điều kiện đầu:    u(x, 0) = sinπx  u(0, t) = 0  u(1, t) = 0          (vd2)  Như vậy với ∆x = xf/M = 1/20 và ∆t = T/N = 1/100 ta có:  ∆t 0.001   r = A 2 = 1. = 0.4               (vd3)  ∆x 0.052 Ta dùng chương trình ctheat.m để tìm nghiệm của (vd1):    clear all, clc  a = 1; %cac thong so cua (vd1)  it0 = inline(ʹsin(pi*x)ʹ,ʹxʹ); %dieu kien dau  bx0 = inline(ʹ0ʹ);   bxf = inline(ʹ0ʹ);%dieu kien bien  xf = 1;   M = 25;   T = 0.1;   N = 100;   [u1, x, t] = fwdeuler(a, xf, T, it0, bx0, bxf, M, N);  figure(1), clf, mesh(t, x, u1)  [u2, x, t] = backeuler(a, xf, T, it0, bx0, bxf, M, N);   figure(2), clf, mesh(t, x, u2)  [u3, x, t] = cranknicholson(a, xf, T , it0, bx0, bxf, M, N);   413
  12. figure(3), clf, mesh(t, x, u3)     4. PDE parabolic 2 chiều: Ta xét bài toán phương trình vi phân đạo hàm riêng  parabolic hai chiều mô tả sự phân bố nhiệt độ u(x, y, t):  ⎡ ∂ 2 u(x,y,t) ∂ 2 u(x,y,t) ⎤ ∂u(x, y,t)   A⎢ + ⎥=           (19)  ⎣ ∂x 2 ∂y 2 ⎦ ∂t Để phương trình có thể giải được ta cần cho điều kiện biên:    u(x0 ,y,t) = bx0 (y,t)   u(xf ,y,t) = b xf (y,t)     u(x, y 0 ,t) = b y0 (x,t)   u(x,y f ,t) = b yf (x,t)   và điều kiện đầu u(x, y, 0) = i0(x, y)  Ta thay đạo hàm bậc 1 theo t ở vế phải bằng sai phân 3 điểm tại điểm giữa  (tk+1 + tk)/2 như phương pháp Crank ‐ Nicholson. Ta cũng thay thế một trong  các đạo hàm bậc hai uxx và uyy bằng xấp xỉ 3 điểm tại thời điểm tk và đạo hàm  kia tại tk+1 và có:  ⎛ u ik,j+1 − 2u ik,j + u ik,j−1 u ik,j+1 − 2u ik,j + u ik,j−1 ⎞ u ik,j+1 − u ik,j   A⎜ − ⎟=      (20)  ⎜ ∆x 2 ∆x 2 ⎟ ∆t ⎝ ⎠ Ta viết phương trình tại thời điểm tiếp theo tk+1:  ⎛ u ik,j++11 − 2u ik,j+1 + u ik,j+−11 u ik,j+1 − 2u ik,j + u ik,j−1 ⎞ u ik,j+ 2 − u ik,j+1   A⎜ − ⎟=     (21)  ⎜ ∆x 2 ∆x 2 ⎟ ∆t ⎝ ⎠ Công  thức  này,  được  Peaceman  và  Rachford  đưa  ra,  là  phương  pháp  ẩn  và  tạo nên hệ phương trình:    ( ) ( ) −ry u ik−+1,j + u ik++1,j + (1 + 2ry )u ik,j+1 = rx u ik,j−1 − u ik,j+1 + (1 − 2rx )u ik,j     1 1 (22a)  với 0 ≤ j ≤ Mx ‐ 1    ( ) ( ) −rx u ik,j+−21 + u ik,j++21 + (1 + 2rx )u ik,j+ 2 = ry u ik−+1,j − u ik++1,j + (1 − 2ry )u ik,j+1   1 1 (22b)  với 0 ≤ i ≤ My ‐ 1  ∆t ∆t và:  rx = A 2     ry = A   ∆x ∆y 2 x − x0 y − y0 T   ∆x = f   ∆y = f   ∆t =   Mx My N Ta xây dựng hàm heat2D() để thực hiện thuật toán này:    function [u, x, y, t] = heat2D(a, D, T, ixy0, bxyt, Mx, My, N)  % Giai au_t = c(u_xx + u_yy) voi D(1) 
  13. % Dieu kien dau: u(x, y, 0) = ixy0(x, y)  % Dieu kien bien: u(x, y, t) = bxyt(x, y, t) voi (x, y)cB  % Mx/My ‐ cac doan co doc theo truc x/y  % N ‐ cac khoang thoi gian  dx = (D(2) ‐ D(1))/Mx;   x = D(1) + [0:Mx]*dx;  dy = (D(4) ‐ D(3))/My;   y = D(3) + [0:My]ʹ*dy;  dt = T/N;   t = [0:N]*dt;  %Khoi gan  for j = 1:Mx + 1      for i = 1:My + 1          u(i, j) = ixy0(x(j), y(i));      end  end  rx = a*dt/(dx*dx);   rx1 = 1 + 2*rx;   rx2 = 1 ‐ 2*rx;  ry = a*dt/(dy*dy);   ry1 = 1 + 2*ry;   ry2 = 1 ‐ 2*ry;  for j = 1:Mx ‐ 1 %Pt.(22a)      Ay(j, j) = ry1;      if j > 1          Ay(j ‐ 1, j) = ‐ry;           Ay(j, j‐1) = ‐ry;       end  end  for i = 1:My ‐ 1 %Pt.(22b)      Ax(i,i) = rx1;      if i > 1          Ax(i ‐ 1, i) = ‐rx;          Ax(i, i ‐ 1) = ‐rx;       end  end  415
  14. for k = 1:N      u_1 = u;       t = k*dt;      for i = 1:My + 1 %Dieu kien bien          u(i, 1) = feval(bxyt, x(1), y(i), t);          u(i, Mx+1) = feval(bxyt, x(Mx+1), y(i), t);      end      for j = 1:Mx + 1          u(1, j) = feval(bxyt, x(j), y(1), t);          u(My+1, j) = feval(bxyt, x(j), y(My + 1), t);      end      if mod(k, 2) == 0          for i = 2:My              jj = 2:Mx;              bx = [ry*u(i, 1) zeros(1, Mx ‐ 3) ry*u(i, My + 1)] ...              +rx*(u_1(i‐1,jj)+ u_1(i + 1,jj)) + rx2*u_1(i,jj);              u(i, jj) = trid(Ay, bxʹ)ʹ; %Pt.(22a)          end      else          for j = 2:Mx              ii = 2:My;              by = [rx*u(1, j); zeros(My‐3,1); rx*u(Mx + 1,j)] ...              + ry*(u_1(ii, j‐1) + u_1(ii, j + 1)) + ry2*u_1(ii, j);              u(ii, j) = trid(Ax, by); %Pt.(22b)          end      end  end    Ta xét phương trình:  −4 ⎡ ∂ u(x,y,t) ∂ 2 u(x,y,t) ⎤ ∂u(x,y,t) 2   10 ⎢ + ⎥=         (vd1)  ⎣ ∂x 2 ∂y 2 ⎦ ∂t trong miền: 0 ≤ x ≤ 4, 0 ≤ y ≤ 4 và trong khoảng thơig gian 0 ≤ t ≤ 5000  Điều kiện đầu:    u(x, y, 0) = 0                  (vd2a)  và điều kiện biên:  eycosx ‐ excosy tại x = 0, x = 4; y = 0, y = 4            (vd2b)  416
  15. Chương trình chương trình ctheat2D.m dùng để giải phương trình là:     clear, clc, clf  a = 1e‐4;  it0 = inline(ʹ0ʹ,ʹxʹ,ʹyʹ); %(vd2a)  bxyt = inline(ʹexp(y)*cos(x)‐exp(x)*cos(y)ʹ,ʹxʹ,ʹyʹ,ʹtʹ); %(vd.2b)  D = [0 4 0 4];   T = 5000;   Mx = 40;   My = 40;   N = 50;  [u, x, y, t] = heat2D(a, D, T, it0, bxyt, Mx, My, N);  mesh(x, y, u)      §4. PHƯƠNG TRÌNH HYPERBOLIC  1.  Dạng  phương  trình:  Phương  trình  truyền  sóng  một  chiều  là  PDE  dạng  hyperbolic:  ∂ 2 u(x,t) ∂ 2 u(x,t)   A =                 (1)  ∂x 2 ∂t 2   0 ≤ x ≤ xf, 0 ≤ t ≤ T  Điều kiện biên:    u(0, t) = b0(t),  u(xf ,t) = b xf (t)   và điều biên:  ∂u   u(x, 0) = i0(x),  = i′ (x)   ∂t t =0 0 phải được cho trước để phương trình có thể giải được  2.  Phương  pháp  sai  phân  tường  minh:  Tương  tự  như  khi  giải  PDE  dạng  parabolic, ta thay đạo hàm bậc hai ở hai vế của (1) bằng sai phân 3 điểm:  u k − 2u ik + u ik−1 u ik +1 − 2u ik + u ik −1     A i +1 =           (2)  ∆x 2 ∆t 2 x T   ∆x = f   ∆t =   M N và có được phương pháp sai phân tường minh:    u ik +1 = r ( u ik+1 + u ik−1 ) + 2(1 − r)u ik − u ik−1           (3)  ∆t 2 với:  r = A   ∆x 2 417
  16. Vì  u i−1 = u(xi , −∆t)  không cho trước nên ta không thể dùng trực tiếp  u1  từ (3)  i với k = 0:    u1 = r ( u i0+1 + u i0−1 ) + 2(1 − r)u i0 − u i−1     i         (4)  Như vậy, ta xấp xỉ điều kiện đầu về đạo hàm bằng sai phân:  u1 − u i−1   i = i′ (xi )                   (5)  2 ∆t 0 và rút ra  u i−1  để đưa vào (3):    u1 = r ( u 0+1 + u 0−1 ) + 2(1 − r)u i0 − ⎡ u1 − 2i′ (xi )∆t ⎤   i i i ⎣ i 0 ⎦ u1 = r ( u 0+1 + u 0−1 ) + (1 − r)u 0 + i′ (xi )∆t     1     i i i i 0       (6)  2 Ta dùng (6) cùng với điều kiện đầu để có  u1  và rồi thay vào (3). Chú ý là:  i   • r ≤ 1 để bảo đảm ổn định  • độ chính xác của nghiêm tăng khi r tăng để cho ∆x giảm  Hợp lí nhất là lấy r = 1. Điều kiện ổn định có thể nhận được bằng cách thay (4)  vào (3):  π 1 λ = 2rcos + 2(1 − r) −   P λ ⎡ π ⎤ hay:  λ 2 + 2 ⎢r ⎛ 1 − cos ⎞ − 1⎥ λ + 1 = 0   ⎜ ⎟ ⎣ ⎝ P⎠ ⎦ Như vậy:  1 ∆t 2   r≤ r = A 2 ≤ 1  ⎛ cos π ⎞′ ∆x 1−⎜ ⎟ ⎝ P⎠ Ta xây dựng hàm wave() để thực hiện thuật toán trên:      function [u, x, t] = wave(a, xf, T, it0, i1t0, bx0, bxf, M, N)  % giai au_xx = u_tt voi  0
  17. for i = 1:M + 1      u(i,1) = it0(x(i));   end  for k = 1:N + 1      u([1 M + 1],k) = [bx0(t(k)); bxf(t(k))];  end  r = a*(dt/dx)^2;   r1 = r/2;   r2 = 2*(1 ‐ r);  u(2:M, 2) = r1*u(1:M ‐ 1, 1) + (1 ‐ r)*u(2:M, 1) + r1*u(3:M + 1, 1) ...               + dt*i1t0(x(2:M)); %Pt.(6)  for k = 3:N + 1      u(2:M, k) = r*u(1:M ‐ 1, k ‐ 1) + r2*u(2:M, k‐1) + r*u(3:M + 1, k ‐ 1)...                     ‐ u(2:M,k ‐ 2); %Pt.(3)  end    Ta xét phương trình:  ∂ 2 u(x,t) ∂ 2 u(x,t)   =                 (vd1)  ∂x 2 ∂t 2   0 ≤ x ≤ 2, 0 ≤ t ≤ 2  Điều kiện đầu và điều kiện biên:  ∂u   u(x, 0) = x(1 ‐ x)     =0            (vd2a)  ∂t t =0 u(0,t) = 0        u(1, t) = 0             (vd2b)  Ta dùng chương trình ctwave.m để giải phương trình này:    clear all, clc  a = 1;  it0 = inline(ʹx.*(1‐x)ʹ,ʹxʹ);  i1t0 = inline(ʹ0ʹ); %(vd2a)  bx0t = inline(ʹ0ʹ);   bxft = inline(ʹ0ʹ); %(vd2b)  xf = 1;   M = 20;   T = 2;   N = 50;  419
  18. [u,x,t] = wave(a, xf, T, it0, i1t0, bx0t, bxft, M, N);  figure(1), clf  mesh(t, x, u)  figure(2), clf  for n = 1:N %hinh dong      plot(x, u(:, n)), axis([0 xf ‐0.3 0.3]), pause(0.2)  end     4. PDE hyperbolic 2 chiều: Phương trình truyền sóng hai chiều là PDE dạng  hyperbolic:  ⎡ ∂ 2 u(x,y,t) ∂ 2 u(x,y,t) ⎤ ∂ 2 u(x,y,t)   A⎢ + ⎥=          (8)  ⎣ ∂x 2 ∂y 2 ⎦ ∂t 2   0 ≤ x ≤ xf, 0 ≤ y ≤ yf, 0 ≤ t ≤ T  Điều kiện biên:    u(0,y,t) = bx0 (y,t)   u(xf , y,t) = b xf (y,t)     u(x,0,t) = b y0 (x,t)  u(x,y f ,t) = b yf (x,t)   và điều biên:  ∂u(x,y)     u(x,y,0) = i 0 (x,y)  = i′ (x, y)   ∂t 0 t =0 Tương tự như hàm một biến, ta thay đạo hàm bậc 2 bằng xấp xỉ 3 điểm:  ⎛ u ik,j+1 − 2u ik,j + u ik,j−1 u ik+1,j − 2u ik,j + u ik−1,j ⎞ u ik,j+1 − 2u ik,j + u ik,j−1     A⎜ + ⎟=   (9)  ⎜ ∆x 2 ∆y 2 ⎟ ∆t 2 ⎝ ⎠ x y T   ∆x = f   ∆y = f   ∆t =   Mx My N và nhận đi đến phương pháp tường minh:  ( ) ( ) u ik,j+1 = rx u ik,j+1 + u ik,j−1 + 2(1 − r x −ry )u ik,j + ry u ik+1,j + u ik−1,j − u ik,j−1   (10)  với:    ∆t 2 ∆t 2   rx = A    ry = A 2   ∆x 2 ∆y Vì  u i−,j = u(xi ,y i , −∆t)  không cho trước nên ta không thể dùng trực tiếp  u1,j  từ  1 i (10) với k = 0:    i (i i ) ( i ) u1,j = rx u 0,j+1 + u 0,j−1 + 2(1 − rx − ry )u i0,j + ry u i0+1,j + u 0−1,j − u i−,j     1 (11)  Như vậy, ta xấp xỉ điều kiện đầu về đạo hàm bằng sai phân:  420
  19. u1,j − u i−,j 1 = i′ (x j ,y i )   i                 (12)  2 ∆t 0 và rút ra  u i−,j  để đưa vào (11):  1   i 1 ( )i ( i i ) u1,j = ⎡rx u i0,j+1 + u 0,j−1 + ry u 0+1,j + u 0−1,j ⎤ 2⎣ ⎦         (13)  + 2(1 − rx − ry )u i0,j + i′ (x j ,y i )∆t 0 Điều kiện ổn định:  4A∆t 2   r= 2 ≤ 1  ∆x + ∆y 2 Ta xây dựng hàm wave2D() để thực hiện thuật toán trên:    function [u,x,y,t] = wave2D(a, D, T, it0, i1t0, bxyt, Mx, My, N)  % giai a(u_xx + u_yy) = u_tt voi D(1) 
  20. rxy1 = 1‐ rx ‐ ry;   rxy2 = rxy1*2;  u_1 = u;  for k = 0:N      t = k*dt;      for i = 1:My + 1 %dieu kien bien          u(i, [1 Mx + 1]) = [bxyt(x(1), y(i),t)  bxyt(x(Mx + 1), y(i),t)];      end      for j = 1:Mx + 1          u([1 My + 1], j) = [bxyt(x(j),y(1),t); bxyt(x(j),y(My + 1),t)];      end      if k == 0          for i = 2:My              for j = 2:Mx %Pt.(13)                  u(i, j) = 0.5*(rx*(u_1(i, j ‐ 1) + u_1(i, j + 1))...                          + ry*(u_1(i ‐ 1,j)+u_1(i + 1,j))) + rxy1*u(i,j) + dt*ut(i,j);              end          end      else          for i = 2:My              for j = 2:Mx                   u(i, j) = rx*(u_1(i, j ‐ 1)+ u_1(i, j + 1))...                          + ry*(u_1(i ‐ 1, j) + u_1(i + 1, j)) + rxy2*u(i, j) ‐ u_2(i, j);              end          end      end      u_2 = u_1;       u_1 = u;       mesh(x, y, u), axis([0 2 0 2 ‐.1 .1])      pause(0.1);  end    Ta xét phương trình:  1 ⎡ ∂ 2 u(x,y,t) ∂ 2 u(x,y,t) ⎤ ∂ 2 u(x,t)   ⎢ + ⎥=           (vd1)  4⎣ ∂x 2 ∂y 2 ⎦ ∂t 2   0 ≤ x ≤ 2, 0 ≤ y ≤ 2, 0 ≤ t ≤ 2  422

CÓ THỂ BẠN MUỐN DOWNLOAD

Đồng bộ tài khoản