Bài giảng Phương trình vi phân đạo hàm riêng
Các PDE được phân thành 3 loại: PDE elliptic: −< 2B4AC0 PDE parabolic: −= 2B4AC0 PDE hyperbolic: −> 2B4AC0 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
Bạn đang xem trước 20 trang tài liệu Bài giảng Phương trình vi phân đạo hàm riêng, để xem tài liệu hoàn chỉnh bạn click vào nút DOWNLOAD ở trên
403
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:
⎛ ⎞∂ ∂ ∂ ∂ ∂+ + = ⎜ ⎟∂ ∂ ∂ ∂ ∂ ∂⎝ ⎠
2 2 2
2 2
u u u u uA(x,y) B(x,y) C(x,y) f x,y,u, ,
x x y y x y
(1)
với xo ≤ x ≤ xf, yo ≤ y ≤ yf và các điều kiện biên:
=o you(x,y ) b (x) =f yfu(x,y ) b (x)
=o xou(x ,y) b (y) =f xfu(x ,y) b (y) (2)
Các PDE được phân thành 3 loại:
PDE elliptic: − <2B 4AC 0
PDE parabolic: − =2B 4AC 0
PDE hyperbolic: − >2B 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 2
2
2 2
u(x,y) u(x,y)u(x,y) g(x,y) g(x,y)u(x,y) f(x,y)
x y
(1)
trên miền { }= ≤ ≤ ≤ ≤o f o fD (x,y) : x x x ,y y y với điều kiện biên dạng:
= =
= =
o yo f yf
o xo f xf
u(x,y ) b (x) u(x,y ) b (x)
u(x ,y) b (y) u(x ,y) b (y)
(2)
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:
+ −
− +∂ ≅∂ ∆
j i
2
i ,j 1 i ,j i ,j 1
2 2
x ,y
u 2u uu(x,y)
x x
với xj = xo + j∆x, yj = yo + j∆y (3.a)
404
+ −− +∂ ≅∂ ∆
j i
2
i 1,j i ,j i 1,j
2 2
x ,y
u 2u uu(x,y)
y x
với ui,j = u(xj, yi) (3.b)
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:
+ − + −
− + − ++ = =∆ ∆
i ,j 1 i ,j i ,j 1 i 1,j i ,j i 1,j
i ,j i ,j i ,j2 2
u 2u u u 2u u
g u f
x y
(4)
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 { }− − − − −x x y y x1,1 1,2 1,M 1 2,1 2,M 1 M 1,2 M 1,M 1u ,u ,...,u ,u ,...,u ,...,u ,...,u . Để
dễ dàng ta viết lại phương trình và điều kiện biên dưới dạng:
+ − + −= + + + + −i ,j y i ,j 1 i ,j 1 x i 1,j i 1,j xy i ,j i ,j i ,ju r (u u ) r (u u ) r (g u f ) (5a)
= = = =
x yi ,o xo i i ,M xf i o,j yo j M ,j yf j
u b (y ) u b (y ) u b (x ) u b (x ) (5b)
Trong đó:
∆= ∆ + ∆
2
y 2 2
yr
2( x y )
∆= ∆ + ∆
2
x 2 2
xr
2( x y )
∆ ∆= ∆ + ∆
2 2
xy 2 2
x yr
2( x y )
(6)
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:
=
∂ ′=∂ o
o
x
x x
u(x,y) b (y)
x
(7)
Thay đạo hàm bậc 1 ở biên trái (x = xo) bằng xấp xỉ 3 điểm:
−
− ′=∆ o
i ,1 i , 1
x i
u u
b (y )
2 x
− ′≈ − ∆oi , 1 i ,1 x iu u 2b (y ) x i = 1, 2,..., My‐1 (8)
Thay thế ràng buộc này vào (5a) ở các điểm biên ta có:
− + −= + + + + −i ,0 y i ,1 i , 1 x i 1,0 i 1,0 xy i ,0 i ,0 i ,0u r (u u ) r (u u ) r (g u f )
+ −′⎡ ⎤= + − ∆ + + + −⎣ ⎦0y i ,1 i ,1 x i x i 1,0 i 1,0 xy i ,0 i ,0 i ,0r u u 2b (y ) x r (u u ) r (g u f )
+ − ′⎡ ⎤= + + + − − ∆⎣ ⎦0y i ,1 x i 1,0 i 1,0 xy i ,0 i ,0 i ,0 x i2r u r (u u ) r g u f 2b (y ) 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:
+ − ⎡ ⎤′= + + + − − ∆⎣ ⎦00,j x 1,j y 0 ,j 1 0,j 1 xy 0,j 0 ,j 0 ,j y ju 2r u r (u u ) r g u f 2b (x ) y (10)
và bổ sung cho góc dưới trái(xo, yo):
405
′′⎡ ⎤= + + − − +⎢ ⎥∆ ∆⎣ ⎦
00 y 0x 0
0,0 y 0,1 x 1,0 xy 0,0 0,0 0 ,0
b (x )b (y )
u 2(r u r u ) r g u f 2 2
x y
(11)
Đ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 <= x <= xf, y0 <= y <= yf}
% voi dieu kien bien:
% u(x0,y) = bx0(y), u(xf,y) = bxf(y)
% u(x,y0) = by0(x), u(x,yf) = byf(x)
% Mx ‐ so doan con tren truc x
% My ‐ so doan con tren truc y
% tol : sai so cho phep
% maxiter: so lan lap
x0 = D(1);
xf = D(2);
y0 = D(3);
yf = D(4);
dx = (xf ‐ x0)/Mx;
x = x0 + [0:Mx]*dx;
dy = (yf ‐ y0)/My;
y = y0 + [0:My]ʹ*dy;
Mx1 = Mx + 1;
My1 = My + 1;
%dieu kien bien
for m = 1:My1
u(m, [1 Mx1]) = [bx0(y(m)) bxf(y(m))];
end
for n = 1:Mx1
u([1 My1], n) = [by0(x(n)); byf(x(n))];
end
sumbv = sum(sum([u(2:My, [1 Mx1]) u([1 My1], 2:Mx)ʹ]));
u(2:My, 2:Mx) = sumbv/(2*(Mx + My ‐ 2));
for i = 1:My
406
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))) < tol
break;
end
u0 = u;
end
Ta giải phương trình Laplace:
∂ ∂∇ = + =∂ ∂
2 2
2
2 2
u(x,y) u(x,y)u(x,y) 0
x y
(vd1)
trong miền 0 ≤ x ≤ 4, 0 ≤ y ≤ 4 với điều kiện biên:
= − = −y x 4u(0,y) e cosy u(4,y) e cos4 e cosy (vd2)
= =x 4 xu(x,0) cosx ‐ e u(x,4) e cosx ‐ e cos4 (vd3)
Ta muốn nhận được u(x, y), mô tả phân bố nhiệt độ trên một tấm vuông mỗi
cạnh dài 4 đơn vị. Ta dùng chương trình ctpoisson.m gọi hàm poisson() để
giải bài toán này.
clear all, clc
407
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
2
u(x,t) u(x,t)A
x t
(1)
Để phương trình có thể giải được ta phải cho điều kiện biên u(0, t) = b0(t),
=
ff x
u(x ,t) b (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
∆ = fx x /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:
+
+ −− + −=∆ ∆
k k k k 1 k
i 1 i i 1 i i
2
u 2u u u uA
x t
(2)
408
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:
+ + −
∆= + + − = ∆
k 1 k k k
i i 1 i 1 i 2
tu r(u u ) (1 2r)u r A
x
(3)
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
k k P
iu e (4)
với P là số nguyên khác zero vào phương trình (3) và có:
π π− π⎡ ⎤⎞⎛λ = + + − = − −⎜ ⎟ ⎢ ⎥⎝ ⎠ ⎣ ⎦
j j
P Pr e e (1 2r) 1 2r 1 cos
P
(5)
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à:
∆= ≤∆ 2
t 1r A
x 2
(6)
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 <= x <= xf, 0 <= t <= T
% dieu kien dau: u(x,0) = it0(x)
ieu kien bien: u(0,t) = bx0(t), u(xf,t) = bxf(t)
% M ‐ so doan con theo x
% N ‐ so diem theo t
dx = xf/M;
x = [0:M]ʹ*dx;
dt = T/N;
t = [0:N]*dt;
for i = 1:M + 1
u(i,1) = it0(x(i));
end
for n = 1:N + 1
u([1 M + 1], n) = [bx0(t(n)); bxf(t(n))];
end
r = a*dt/dx/dx
r1 = 1 ‐ 2*r;
for k = 1:N
for i = 2:M
u(i, k+1) = r*(u(i + 1, k) + u(i‐1, k)) + r1*u(i, k); %Pt.(3)
409
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):
−
+ −− + −=∆ ∆
k k k k k 1
i 1 i i 1 i i
2
u 2u u u uA
x t
(7)
−− +
∆− + + − = = ∆
k k k k 1
i 1 i i 1 i 2
tru (1 2r)u ru u r A
x
(8)
Nếu các giá trị k0u và
k
Mu ở 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:
−
−
−
−
− −
−
− −
⎡ ⎤ ⎡ ⎤+⎡ ⎤+ − ⎢ ⎥ ⎢ ⎥⎢ ⎥− + − ⎢ ⎥ ⎢ ⎥⎢ ⎥ ⎢ ⎥ ⎢ ⎥⎢ ⎥− + − ⎢ ⎥ ⎢ ⎥⎢ ⎥ =⎢ ⎥ ⎢ ⎥⎢ ⎥ ⎢ ⎥ ⎢ ⎥⎢ ⎥+ − ⎢ ⎥ ⎢ ⎥⎢ ⎥ ⎢ ⎥ ⎢ ⎥⎢ ⎥− + +⎣ ⎦ ⎣ ⎦ ⎣ ⎦
M M M M M M M M
L
L
k k 1 k
1 1 0
k k 1
2 2
k k 1
3 3
k k 1
M 2 M 2
k k 1 k
M 1 M 1 M
u u ru1 2r r 0 0 0 0
u ur 1 2r r 0 0 0
0 r 1 2r r 0 0 u u
0 0 0 1 2r r u u
0 0 0 r 1 2r u u ru
(9)
Điều kiện biên Neumann
=
∂ ′=∂ 0x 0
u b (t)
x
được đưa vào phương trình bằng cách
xấp xỉ:
−− ′=∆
k k
1 1
0
u u b (k)
2 x
(10)
và ghép nó với phương trình có ẩn k0u :
−−− + + − =k k k k 11 0 1 0ru (1 2r)u ru u (11)
để có được phương trình:
− ′+ − = − ∆k k k 10 1 0 0(1 2r)u 2ru u 2rb (k) x (12)
Kết quả ta có được hệ phương trình:
410
−
−
−
−
−
−
−
⎡ + − ′⎡ − ∆⎡ ⎤⎤⎢ ⎢⎢ ⎥⎥− + −⎢ ⎢⎢ ⎥⎥⎢ ⎢⎢ ⎥⎥− + −⎢ ⎢⎢ ⎥⎥− + =⎢ ⎢⎢ ⎥⎥⎢ ⎢⎢ ⎥⎥−⎢ ⎢⎢ ⎥⎥+ −⎢ ⎢⎢ ⎥⎥⎢ ⎢⎢ ⎥⎥− + +⎦ ⎣ ⎦ ⎣⎣
L
L
L
L
M M M M M M M
L
L
k k
0 0 0
k k 1
1 1
k k 1
2 2
k k 1
3 3
k k 1
0 M 2
k k 1 k
0 M 1 M
1 2r r 0 0 0 0 u u 2rb (k) x
r 1 2r r 0 0 0 u u
0 r 1 2r r 0 0 u u
0 0 r 1 2r 0 0 u u
r 0
0 0 0 0 1 2r r u u
0 0 0 0 r 1 2r u u ru
⎤⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦
(13
Điểu kiện ổn định của nghiệm là:
π π
−− + + − = λ
j j
P P 1re (1 2r) re
hay: λ = π⎡ ⎤+ −⎢ ⎥⎣ ⎦
1
1 2r 1 cos
P
λ ≤ 1 (14)
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 <= x <= xf, 0 <= t <= T
% Dieu kien dau: u(x,0) = it0(x)
% ieu kien bien: u(0,t) = bx0(t), u(xf,t) = bxf(t)
% M ‐ so khoang con tren truc x
% N ‐ so khoang theo t
dx = xf/M;
x = [0:M]ʹ*dx;
dt = T/N;
t = [0:N]*dt;
for i = 1:M + 1
u(i, 1) = it0(x(i));
end
for n = 1:N + 1
u([1 M + 1], n) = [bx0(t(n)); bxf(t(n))];
end
r = a*dt/dx/dx;
r2 = 1 + 2*r;
for i = 1:M ‐ 1
411
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ó:
+ + + +
+ − + −⎛ ⎞− + − + −+ =⎜ ⎟∆ ∆ ∆⎝ ⎠
k 1 k 1 k 1 k k k k 1 k
i 1 i i 1 i 1 i i 1 i i
2 2
A u 2u u u 2u u u u
2 x x t
(15)
và nhận được phương pháp Crank ‐ Nicholson:
+ + ++ − + −
∆− + + − = + − + = ∆
k 1 k 1 k 1 k k k
i 1 i i 1 i 1 i i 1 2
tru (1 2r)u ru ru (1 2r)u ru r A
x
(16)
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:
+
+
+
+
−
+
⎡ ⎤⎡ ⎤+ − ⎢ ⎥⎢ ⎥− + − ⎢ ⎥⎢ ⎥ ⎢ ⎥⎢ ⎥− + − ⎢ ⎥⎢ ⎥ ⎢ ⎥⎢ ⎥ ⎢ ⎥⎢ ⎥+ − ⎢ ⎥⎢ ⎥ ⎢ ⎥⎢ ⎥− +⎣ ⎦ ⎣ ⎦
M M M M M M M
L
L
k 1
1
k 1
2
k 1
3
k 1
M 1
k 1
M
u2(1 r) r 0 0 0 0
ur 2(1 r) r 0 0 0
0 r 2(1 r) r 0 0 u
0 0 0 2(1 r) r u
0 0 0 r 2(1 r) u
−
⎡ ⎤⎡ ⎤− ⎢ ⎥⎢ ⎥− ⎢ ⎥⎢ ⎥ ⎢ ⎥⎢ ⎥− ⎢ ⎥⎢ ⎥= ⎢ ⎥⎢ ⎥ ⎢ ⎥⎢ ⎥− ⎢ ⎥⎢ ⎥ ⎢ ⎥⎢ ⎥−⎣ ⎦ ⎣ ⎦
M M M M M M M
L
L
k
1
k
2
k
3
k
M 1
k
M
u2(1 r) r 0 0 0 0
ur 2(1 r) r 0 0 0
0 r 2(1 r) r 0 0 u
0 0 0 2(1 r) r u
0 0 0 r 2(1 r) u
412
[ ]
+⎡ ⎤+⎢ ⎥⎢ ⎥⎢ ⎥⎢ ⎥+ ⎢ ⎥⎢ ⎥⎢ ⎥′ ′⎢ ⎥+ +⎣ ⎦
M
k 1 k
0 0
M M
r(u u )
0
0
0
2r b (k 1) b (k)
(17)
Điều kiện ổn định được xác định bằng:
π π⎡ ⎤ ⎡ ⎤⎛ ⎞ ⎛ ⎞λ + − = − −⎜ ⎟ ⎜ ⎟⎢ ⎥ ⎢ ⎥⎝ ⎠ ⎝ ⎠⎣ ⎦ ⎣ ⎦2 1 r 1 cos 2 1 r 1 cosP P
hay:
π⎡ ⎤− −⎢ ⎥⎣ ⎦λ = λ ≤π⎡ ⎤+ −⎢ ⎥⎣ ⎦
1 r 1 cos
P 1
1 r 1 cos
P
(18)
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 <= x <= xf, 0 <= t <= T
% Dieu kien dau: u(x,0) = it0(x)
% Dieu kien bien: u(0, t) = bx0(t), u(xf, t) = bxf(t)
% M ‐ so khoang con tren truc x
% N ‐ so khoang theo t
dx = xf/M;
x = [0:M]ʹ*dx;
dt = T/N;
t = [0:N]*dt;
for i = 1:M + 1
u(i, 1) = it0(x(i));
end
for n = 1:N + 1
u([1 M + 1], n) = [bx0(t(n)); bxf(t(n))];
end
r = a*dt/dx/dx;
r1 = 2*(1 ‐ r);
r2 = 2*(1 + r);
for i = 1:M ‐ 1
413
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
2
u(x,t) u(x,t) 0 x 1, 0 t 0.1
x t
(vd1)
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ó:
∆= = =∆ 2 2
t 0.001r A 1. 0.4
x 0.05
(vd3)
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);
414
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 2
2 2
u(x,y,t) u(x,y,t) u(x,y,t)A
x y t
(19)
Để phương trình có thể giải được ta cần cho điều kiện biên:
=
00 x
u(x ,y,t) b (y,t) =
ff x
u(x ,y,t) b (y,t)
=
00 y
u(x,y ,t) b (x,t) =
ff y
u(x,y ,t) b (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ó:
+
+ − + −⎛ ⎞− + − + −− =⎜ ⎟⎜ ⎟∆ ∆ ∆⎝ ⎠
k k k k k k k 1 k
i ,j 1 i ,j i ,j 1 i ,j 1 i ,j i ,j 1 i ,j i ,j
2 2
u 2u u u 2u u u u
A
x x t
(20)
Ta viết phương trình tại thời điểm tiếp theo tk+1:
+ + + + +
+ − + −⎛ ⎞− + − + −− =⎜ ⎟⎜ ⎟∆ ∆ ∆⎝ ⎠
k 1 k 1 k 1 k k k k 2 k 1
i ,j 1 i ,j i ,j 1 i ,j 1 i ,j i ,j 1 i ,j i ,j
2 2
u 2u u u 2u u u u
A
x x t
(21)
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:
( ) ( )+ + +− + − +− + + + = − + −k 1 k 1 k 1 k k ky i 1,j i 1,j y i ,j x i ,j 1 i ,j 1 x i ,jr u u (1 2r )u r u u (1 2r )u (22a)
với 0 ≤ j ≤ Mx ‐ 1
( ) ( )+ + + + + +− + − +− + + + = − + −k 2 k 2 k 2 k 1 k 1 k 1x i ,j 1 i ,j 1 x i ,j y i 1,j i 1,j y i ,jr u u (1 2r )u r u u (1 2r )u (22b)
với 0 ≤ i ≤ My ‐ 1
và: ∆= ∆x 2
tr A
x
∆= ∆y 2
tr A
y
−∆ f 0
x
x xx=
M
−∆ f 0
y
y yy=
M
∆ = Tt
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) <= x <= D(2), D(3) <= y <= D(4), 0 <= t
%<= T
415
% 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
416
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:
−
⎡ ⎤∂ ∂ ∂+ =⎢ ⎥∂ ∂ ∂⎣ ⎦
2 2
4
2 2
u(x,y,t) u(x,y,t) u(x,y,t)10
x y t
(vd1)
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)
417
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 2
2 2
u(x,t) u(x,t)A
x t
(1)
0 ≤ x ≤ xf, 0 ≤ t ≤ T
Điều kiện biên:
u(0, t) = b0(t), =
ff x
u(x ,t) b (t)
và điều biên:
u(x, 0) = i0(x),
=
∂ ′=∂ 0t 0
u i (x)
t
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:
+ −
+ −− + − +=∆ ∆
k k k k 1 k k 1
i 1 i i 1 i i i
2 2
u 2u u u 2u uA
x t
(2)
∆ = fxx
M
∆ = Tt
N
và có được phương pháp sai phân tường minh:
( )+ −+ −= + + − −k 1 k k k k 1i i 1 i 1 i iu r u u 2(1 r)u u (3)
với: ∆= ∆
2
2
tr A
x
418
Vì − = −∆1i iu u(x , t) không cho trước nên ta không thể dùng trực tiếp 1iu từ (3)
với k = 0:
( ) −+ −= + + − −1 0 0 0 1i i 1 i 1 i iu r u u 2(1 r)u u (4)
Như vậy, ta xấp xỉ điều kiện đầu về đạo hàm bằng sai phân:
−− ′=∆
1 1
i i
0 i
u u i (x
