Giáo trình Matlab - Chương 8: Phương trình vi phân đạo hàm riêng (Phần 1)

ành: 
0ucu2 =∆−ω−  
160
hay phương trình Helmholz: 
0uku 2 =−∆−  
Trong đó k  là số sóng,  liên quan đến  tần số góc ω,  tần số  f và độ dài sóng λ 
bằng: 
λ
π=π=ω= 2
c
f2
c
k  
Bây giờ  ta phải mô  tả  điều kiện biên. Coi  sóng  tới  là một  sóng phẳng  theo 
hướng: 
))asin(),a(cos(a =r : 
ti)tx.ak(i e)y,x(ve)t,y,x(v ω−ω− == rr  
Trong đó :   )y,x(ve x.aik =rr
u là tổng của v và sóng phản xạ r: 
u = v + r 
Điều kiện biên đối với biên đối tượng là u = 0, nghĩa là r = ‐ v(x,y)(với sóng âm, 
v là nhiễu loạn áp suất do đó điều kiện biên thích hợp phải là ∂u/∂n=0). Sóng 
phản xạ đi  ra khỏi  đối  tượng. Điều kiện biên bên ngoài phải  chọn  sao  sóng 
không bị phản xạ. Điều kiện như vậy thường gọi là không phản xạ và ta dùng 
điều kiện bức xạ Sommerfeld. Khi đạt đến vô cùng, r thoả mãn xấp xỉ phương 
trình sóng một chiều: 
0rc
t
r =∇ξ+∂
∂ r  
cho phép sóng chỉ chuyển động theo hướng x(x là khoảng cách bức xạ từ vật 
thể). Với nghiệm điều hoà theo t,đ iều này trở thành điều kiện biên Neumann 
tổng quát hoá: 
ikrr =∇ξr  
Vì lí do đơn giản hoá, coi pháp tuyến bên ngoài của vùng khảo sát xấp xỉ 
hướng bên ngoài ξ. 
  Bây giờ ta dùng pdetool để giải bài toán. Sử dụng mode Generic Scalar 
bằng cách vào menu Option|Application và đánh dấu vào Generic Scalar , bắt 
đầu bằng cách vẽ vùng 2‐D của bài toán.  
Đối tượng bị chiếu sáng là một hình vuông R1 với cạnh 0.1 đơn vị và tâm 
ở [0.8 0.5] và quay 450 (vào Draw|Rotate ) và coi vùng tính là C1 là hình tròn 
có bán kính 0.45 và tâm cùng một chỗ với hình vuông. Mô hình CGS là C1‐R1.  
Với biên ngoài, ta dùng điều kiện biên Neumann với q = ‐jk. Hệ số sóng 
k = 60 tương ứng với bước sóng 0.1 đơn vị. Như vậy ta cần nhập giá trị q = ‐60j 
và g = 0 bằng cách vào menu Boundary| Specify Boundary Condition...  
Đối với biên hình vuông ta dùng điều kiện biên Dirichlet: 
161
x.aike)y,x(vr
rr−=−=  
 Trong bài  toán này,  sóng  tới phải  đi qua  đoạn  ‐x. Như vậy  điều kiện 
biên có dạng đơn giản là: 
ikxe)y,x(vr −=−=  
Vào menu Boundary| Specify Boundary Condition... và nhập điều kiện 
biên Dirichlet h = 1, r =  ‐exp(‐i*60*x). Vào menu PDE|PDE Specification... và 
nhập các hệ số của phương trình là c = 1, a = ‐3600 và f = 0. Tiếp theo bấm vào 
nút công cụ chia  lưới để chia  lưới cho bài toán. Ta  lưu tất cả kết quả vào file 
ct8_4.m. 
Bây giờ ta có thể giải bài toán và có nghiệm phức.  
Để thấy được sự  lan truyền sóng phản xạ, sau khi giải bài toán ta chọn 
Draw|Export Geometry Description, Set Formular, Label...; Boundary|Export 
Decomposed  Gemetry,  Boundary  Condʹs...;  PDE|Export  Coeficients...; 
Solve|Export  Solution...    và  chạy  các  lệnh  MATLAB  sau(lưu  trong  file 
ct8_5.m): 
m = 10; % so khung hinh 
h = newplot;  
hf = get(h,ʹParentʹ);  
set(hf,ʹRendererʹ,ʹzbufferʹ) 
axis tight; 
set(gca,ʹDataAspectRatioʹ,[1 1 1]);  
axis off; 
M = moviein(m,hf); 
maxu = max(abs(u)); 
for j = 1:m 
  uu = real(exp(‐j*2*pi/m*sqrt(‐1))*u); 
  fprintf(ʹ%d ʹ,j); 
  pdeplot(p,e,t,ʹxydataʹ,uu,ʹcolorbarʹ,ʹoffʹ,ʹmeshʹ,ʹoffʹ); 
  caxis([‐maxu maxu]); 
  axis tight, set(gca,ʹDataAspectRatioʹ,[1 1 1]);  
  axis off; 
  M(:,j) = getframe(hf); 
  if j = = m 
    fprintf(ʹdone\nʹ); 
  end 
end 
movie(hf,M,50); 
162
c. Bài toán mặt cực tiểu: Trong nhiều bài toán hệ số c, a và f không chỉ 
phụ thuộc vào x và y mà còn vào u. Ta khảo sát phương trình: 
  0u
|u|1
1.
2
=⎟⎟⎠
⎞
⎜⎜⎝
⎛
∇
∇+
∇−  
trong hình tròn đơn vị  ( ){ }1yx|y,x 22 ≤+=Ω  với u = x2 trên ∂Ω. Đây là bài toán 
phi  tuyến  tính và không  thể giải với  trình giải solver. Do đó  trình pdenonlin 
được dùng. 
Vùng tìm nghiệm là hình tròn đơn vị. Vẽ hình tròn này và xác định điều 
kiện biên bằng cách vào menu Boundary|Boundary Mode. Dùng Select All của 
menu Edit để chọn tất cả các biên. Sau đó nhấp đúp vào biên để mở hộp thoại 
Boundary Condition. Điều kiện biên u = x^2 được nhập vào bằng cách đánh 
x.^2 vào mục  của  r. Tiếp  theo mở hộp  thoại PDE Specification  để xác  định 
PDE. Đây là phương trình elliptic với 
2|u|1
1c ∇+= , a = 0 và f = 0. Hệ số c được 
nhập vào bằng cách đánh c =1./sqrt(1+ux.^2+uy.^2). Khởi gán lưới và làm tinh 
lại một lần. Trước khi giải phương trình chọn Parameters . . . từ menu Solve và 
dùng tuỳ chọn Use nonlinear và đặt sai số 0.001. Nhấn nút = để giải phương 
trình. Dùng hộp thoại Plot Selection để vẽ nghiệm dạng 3‐D. 
d. Chia vùng: PDE Toolbox được thiết kế để làm việc với sự phân vùng 
cấp 1. Nếu vùng Ω phức tạp, thường nên phân nó thành các vùng có cấu trúc 
đơn  giản  hơn  đơn  giản. Các  cấu  trúc  như  vậy  thường  được  thực  hiện  bởi 
pdetool. 
Giả sử Ω là một tập hợp các vùng con Ω1, Ω2,..., Ωn.Vậy ta có thể đánh số 
lại  các nút  của  lưới  trên Ω  sao  cho  chỉ  số  của  các nút  ở  các vùng  con  được 
nhóm  lại với nhau  trong khi  tất cả chỉ số của các nút chung của 2 hay nhiều 
vùng con như cũ. Do ma trận K có các số hạng khác 0 chỉ ở các hàng và các cột 
chỉ các nút cạnh nhau, nên nó được phân thành: 
⎟⎟
⎟⎟
⎟⎟
⎠
⎞
⎜⎜
⎜⎜
⎜⎜
⎝
⎛
⋅⋅⋅
⋅⋅⋅
=
CBBB
BK00
B0K0
B00K
K
n21
T
nn
T
22
T
11
L
L
MMOMM   
Trong khi vế phải là: 
163
⎟⎟
⎟⎟
⎟⎟
⎠
⎞
⎜⎜
⎜⎜
⎜⎜
⎝
⎛
=
c
n
2
1
f
f
f
f
f M  
  Thường trình assempde của PDE Toolbox có thể gộp các ma trận Kj, Bj, fj 
và C một cách riêng rẽ nên ta có thể kiểm soát được việc lươ và xử lí các ma 
trận này. Hơn nữa cấu trúc của hệ thống tuyến tính Ku = F được đơn giản hoá 
bằng cách chia K thành từng ma trận riêng phần như trên. 
  Bây giờ ta khảo sát một màng hình L. Ta có thể vẽ màng này bằng lệnh:
  pdegplot(ʹlshapegʹ) 
Chú ý các biên giữa các vùng con. Có 3 vùng con vì miền đang xét có dạng L. 
Như vậy công thức ma trận với n = 3 từ trên có thể dùng. Bây giờ ta tạo lưới : 
[p,e,t] = initmesh(ʹlshapegʹ); 
[p,e,t] = refinemesh(ʹlshapegʹ,p,e,t); 
[p,e,t] = refinemesh(ʹlshapegʹ,p,e,t); 
Với trường hợp này với n = 3 ta có: 
⎟⎟
⎟⎟
⎟
⎠
⎞
⎜⎜
⎜⎜
⎜
⎝
⎛
=
⎟⎟
⎟⎟
⎟
⎠
⎞
⎜⎜
⎜⎜
⎜
⎝
⎛
⎟⎟
⎟⎟
⎟
⎠
⎞
⎜⎜
⎜⎜
⎜
⎝
⎛
c
3
2
1
3
2
1
321
T
33
T
22
T
11
f
f
f
f
c
u
u
u
CBBB
BK00
B0K0
B00K
Và nghiệm xác định bằng cách loại trừ khối: 
L
)uBf(Ku
fKBfKBfKBfu)BKBBKBBKBC(
c
T
11
1
11
3
1
332
1
221
1
11cc
T
3
1
33
T
2
1
22
T
1
1
11
−=
−−−=−−−
−
−−−−−−
Khi  giải  bài  toán,  ta dùng phương pháp phân  tích Choleski. Cac  lệnh 
MATLAB như sau( lưu trong ct8_7.m): 
echo on 
clc 
%       Gia pt Poissonʹs  
%        ‐div(grad(u))=1 
%       tren mang L voi u=0 tren bien. 
%       Bai toan duoc gia bang cach chia vung. 
pause % Nhan phim bat ki de tiep tuc. 
164
clc 
%       Dinh nghia bai toan 
g = ʹlshapegʹ; % mang dang L 
b = ʹlshapebʹ; % 0 tren bien 
c = 1; 
a = 0; 
f = 1; 
time = [];  
[p,e,t] = initmesh(g); 
[p,e,t] = refinemesh(g,p,e,t); 
[p,e,t] = refinemesh(g,p,e,t); 
pause % Nhan phim bat ki de tiep tuc. 
clc 
np = size(p,2); 
%       Truoc het tim cac diem chung 
cp = pdesdp(p,e,t); 
%       Dinh vi khong gian 
nc = length(cp); 
C = zeros(nc,nc);  
FC = zeros(nc,1); 
pause % Nhan phim bat ki de tiep tuc. 
%       Ket hop vung 1 va cap nhat 
[i1,c1] = pdesdp(p,e,t,1); 
ic1 = pdesubix(cp,c1); 
[K,F] = assempde(b,p,e,t,c,a,f,time,1); 
K1 = K(i1,i1); 
d = symmmd(K1); 
i1 = i1(d); 
K1 = chol(K1(d,d)); 
B1 = K(c1,i1);  
a1 = B1/K1; 
C(ic1,ic1) = C(ic1,ic1) + K(c1,c1) ‐ a1*a1ʹ; 
165
f1 = F(i1); 
e1 = K1ʹ\f1; 
FC(ic1) = FC(ic1) + F(c1) ‐ a1*e1; 
pause % Nhan phim bat ki de tiep tuc. 
%       Ket hop vung 2 va cap nhat 
[i2,c2] = pdesdp(p,e,t,2); 
ic2 = pdesubix(cp,c2); 
[K,F] = assempde(b,p,e,t,c,a,f,time,2); 
K2 = K(i2,i2); 
d = symmmd(K2); 
i2 = i2(d); 
K2 = chol(K2(d,d)); 
B2 = K(c2,i2); 
a2 = B2/K2; 
C(ic2,ic2) = C(ic2,ic2) + K(c2,c2) ‐ a2*a2ʹ; 
f2 = F(i2); 
e2 = K2ʹ\f2; 
FC(ic2) = FC(ic2) + F(c2) ‐ a2*e2; 
pause % Nhan phim bat ki de tiep tuc. 
%       Ket hop vung 3 va cap nhat 
[i3,c3] = pdesdp(p,e,t,3); 
ic3 = pdesubix(cp,c3); 
[K,F] = assempde(b,p,e,t,c,a,f,time,3); 
K3 = K(i3,i3); 
d = symmmd(K3); 
i3 = i3(d); 
K3 = chol(K3(d,d)); 
B3 = K(c3,i3); 
a3 = B3/K3; 
C(ic3,ic3) = C(ic3,ic3) + K(c3,c3) ‐ a3*a3ʹ; 
f3 = F(i3); 
e3 = K3ʹ\f3; 
FC(ic3) = FC(ic3 )+ F(c3)‐a3*e3; 
pause % Nhan phim bat ki de tiep tuc. 
166
% Giai 
u = zeros(np,1); 
u(cp) = C\FC; % cac diem chung 
u(i1) = K1\(e1‐a1ʹ*u(c1)); % cac diem tren vung 1 
u(i2) = K2\(e2‐a2ʹ*u(c2)); % cac diem tren vung 2 
u(i3) = K3\(e3‐a3ʹ*u(c3)); % cac diem tren vung 3 
% Ve 
pdesurf(p,t,u) 
pause % Nhan phim bat ki de tiep tuc. 
echo off 
Các lệnh giải bằng GUI lưu trong ct8_8.m 
2. Các ví dụ về bài toán parabolic: 
a.  Phương  trình  truyền  nhiệt:  khối  kim  loại  bị  đốt  nóng:  Bài  toán 
parabolic nói chung là bài toán về phương trình truyền nhiệt dạng: 
0u
t
ud =∆−∂
∂  
Phương trình truyền nhiệt mô tả sự lan truyền nhiệt trong các vật thể. Ví 
dụ đầu tiên mà ta xem xét là một thanh kim loại có một lỗ hình chữ nhật bị đốt 
nóng. Tại đầu bên trái của thanh nhiệt độ là 1000C. Tại đầu bên phải nhiệt lan 
truyền từ thanh vào môi trường không khí chung quanh với tốc độ hằng, ví dụ 
là 10. Các biên khác cách nhiệt. Do đó ta có các điều kiện biên sau: 
  • u = 100 ở bên trái (điều kiện biên Dirichlet) 
  •  10
n
u −=∂
∂  ở bên phải (điều kiện biên Neumann) 
  •  0
n
u =∂
∂  trên các biên khác(điều kiện biên Neumann) 
Với phương trình truyền nhiệt ta cần điều kiện đầu: nhiệt độ của thanh 
kim  loại  tại  t0. Trong bài  toán này nhiệt  độ ban  đầu  của  thanh  là 00C. Cuối 
cùng ta xác định thời điểm bắt đầu truyền nhiệt là t = 0 và ta muốn tính toán 
quá trình truyền nhiệt trong suốt 5s. 
  Sau khi đã khởi động pdetool và chọn Generic Scalar ta vẽ mô hình CSG  
dễ dàng. Ta vẽ hình chữ nhật R1: [‐0.5 ‐0.8 1 1.6]. Ta vẽ tiếp hình chữ nhật khác 
(lỗ) R2: [‐0.05 ‐0.4  0.1  0.8]. Mô hình này là R1 ‐ R2. 
  Bấm vào nút ∂Ω để xác định biên để mô tả  biên và  điều kiện biên. Dùng 
167