Xem mẫu
- 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 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
- ⎡ 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
- 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)))
- 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
- 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
- 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
-
⎡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
- 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
- ⎡ 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
- 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
- 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)
- % 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
- 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
- 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
- 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
- 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
- [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
- 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)
- 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
nguon tai.lieu . vn