Giáo trình Matlab cơ bản - Chương 9: Phương trình vi phân đạo hàm riêng

1 11 12;1 11 19;10 11 19;4 5 19;5 7 19; 5 6 7;1 2 15; 2 3 15; 
       3 15 17;3 4 17;4 17 19;13 17 19;1 13 19;1 13 15;7 8 22;8 9 22; 
       9 22 24;9 10 24; 10 19 24; 19 20 24;7 19 20; 7 20 22;13 14 18; 
       14 15 16;16 17 18;20 21 25;21 22 23;23 24 25;14 26 28; 
       16 26 27;18 27 28; 21 29 31;23 29 30;25 30 31; 
       26 27 28; 29 30 31]; %mien con tam giac 
fexemp = ʹ(norm([x y] + [0.5 0.5])<0.01) ‐ (norm([x y] ‐ [0.5 0.5]) < 0.01)ʹ; 
f = inline(fexemp,ʹxʹ,ʹyʹ); %(Pt.2) 
g = inline(ʹ0ʹ,ʹxʹ,ʹyʹ); 
Nn = size(N, 1); %tong so nut 
Ni = Nn ‐ Nb; %so nut ben trong 
c = zeros(1, Nn); %gia tri tren bien 
p = fembasisftn(N, S); 
[U, c] = femcoef(f, g, p, c, N, S, Ni); 
%do thi luoi tam giac 
figure(1); 
clf; 
trimesh(S, N(:, 1), N(:, 2), c); 
%do thi luoi chu nhat 
Ns = size(S, 1); %tong so mien con tam giac 
x0 = ‐1;  
xf = 1;  
y0 = ‐1;  
yf = 1; 
 432
Mx = 16;  
dx = (xf ‐ x0)/Mx;  
xi = x0 + [0:Mx]*dx; 
My = 16;  
dy = (yf ‐ y0)/My;  
yi = y0 + [0:My]*dy; 
for i = 1:length(xi) 
    for j = 1:length(yi) 
        for s = 1:Ns  
            if inpolygon(xi(i), yi(j), N(S(s,:), 1), N(S(s,:),2)) > 0 
                Z(i, j) = U(s,:)*[1 xi(i) yi(j)]ʹ; %Pt.(4.5b) 
                break; 
            end 
        end 
    end 
end 
figure(2); 
clf; 
mesh(xi, yi, Z) 
%de so sanh 
bx0 = inline(ʹ0ʹ);  
bxf = inline(ʹ0ʹ); 
by0 = inline(ʹ0ʹ);  
byf = inline(ʹ0ʹ); 
D = [x0 xf y0 yf]; 
[U, x, y] = poisson(f, g, bx0, bxf, by0, byf, D, Mx, My, 1e‐6, 50); 
figure(3) 
clf; 
mesh(x,y,U) 
§6. GUI CỦA MATLAB ĐỂ GIẢI PDE 
1. Các phương trình có thể giải được bằng PDETOOL: Công cụ PDETOOL 
của MATLAB có thể dùng để giải các loại phương trình sau: 
  a. Phương trình elliptic: Ta sẽ giải phương trình elliptic 
  −∇ ∇ + =(c u) au f                   (1) 
với điều kiện bên: 
 433
=
∇r
hu r Dirichlet
nc u+qu=g Neumann
              (2) 
trên biên ∂Ω, trong đó  rn  là vec tơ pháp tuyến. 
  Trong trường hợp u là đại lượng vô hướng, phương trình (1) trở thành: 
⎡ ⎤∂ ∂− + + =⎢ ⎥∂ ∂⎣ ⎦
2 2
2 2
u(x,y) u(x,y)c au(x,y) f(x,y)
x y
        (3) 
và nếu điều kiện biên đối với phân biên bên trái là điều kiện biên Neumann 
dạng 
=
∂ ′=∂ 0
0
x
x x
u(x,y) b (y)
x
 thì (2) có thể viết thành: 
⎡ ⎤∂ ∂− + +⎢ ⎥∂ ∂⎣ ⎦
∂= − + =∂
r r r
x x y
u(x,y) u(x,y)e c e e qu(x,y)
x y
u(x,y)c qu(x,y) g(x,y)
x
          (4) 
vì vec tơ pháp tuyến của biên phải là  =r rxn e  
  b. Phương trình parabolic: Ta xé giải phương trình: 
  ∂−∇ ∇ + + =∂
u(c u) au d f
t
              (5) 
trên miền Ω và trong khoảng thời gian 0 ≤ t ≤ T, với điều kiện bên giống (2) và 
điều kiện đầu u(t0) 
  c. Phương trình hyperbolic: 
  ∂−∇ ∇ + + =∂
2
2
u(c u) au d f
t
              (6) 
trên miền Ω và trong khoảng thời gian 0 ≤ t ≤ T, với điều kiện bên giống (2) và 
điều kiện đầu u(t0), u’(t0) 
  d. Phương trình giá trị riêng: 
∂−∇ ∇ + = λ ∂
u(c u) au
t
                (7) 
trên miền Ω với một giá trị riêng chưa biết λ và điều kiện biên tương tự (2). 
  Công cụ PDETOOL cũng có thể dùng để giải hệ phương trình dạng: 
−∇ ∇ −∇ ∇ + + =⎧⎨−∇ ∇ −∇ ∇ + + =⎩
11 1 12 2 11 1 12 2 1
21 1 22 2 21 1 22 2 2
(c u ) (c u ) a u a u f
(c u ) (c u ) a u a u f
        (8) 
trên miền Ω với điều kiện biên Dirichlet: 
  ⎡ ⎤ ⎡ ⎤ ⎡ ⎤=⎢ ⎥ ⎢ ⎥ ⎢ ⎥⎣ ⎦ ⎣ ⎦⎣ ⎦
1 111 12
2 221 22
u rh h
u rh h
                (9) 
hay điều kiện Neumann tổng quát: 
 434
∇ + ∇ + + =⎧⎨ ∇ + ∇ + + =⎩
r r
r r11 1 12 2 11 1 12 2 1
21 1 22 2 21 1 22 2 2
n(c u ) n(c u ) q u q u g
n(c u ) n(c u ) q u q u g
        (10) 
hay điều kiện biên hỗn hợp: 
  ⎡ ⎤= ⎢ ⎥⎣ ⎦
11 12
21 22
c c
c
c c
  ⎡ ⎤= ⎢ ⎥⎣ ⎦
11 12
21 22
a a
a
a a
  ⎡ ⎤= ⎢ ⎥⎣ ⎦
1
2
f
f
f
  ⎡ ⎤= ⎢ ⎥⎣ ⎦
1
2
u
u
u
  ⎡ ⎤= ⎢ ⎥⎣ ⎦
11 12
21 22
h h
h
h h
  ⎡ ⎤= ⎢ ⎥⎣ ⎦
1
2
r
r
r
  ⎡ ⎤= ⎢ ⎥⎣ ⎦
11 12
21 22
q q
q
q q
  ⎡ ⎤= ⎢ ⎥⎣ ⎦
1
2
g
g
g
2. Sử dụng PDETOOL: PDETOOL giải phương trình vi phân đạo hàm riêng 
bằng cách dùng phương pháp FEM. Để giải phương  trình  ta  theo các bước 
sau: 
) Nhập  lệnh pdetool vào  cửa  sổ  lệnh MATLAB. Cửa  sổ PDE  toolbox 
xuất hiện. Ta có thể bật/tắt tuỳ chọn Grid bằng cách bấm vào Grid trên 
menu Option. Ta cũng có thể hiệu chỉnh phạm vi trục x và y bằng cách 
chọn Axes Limit trong nemu Option 
Nếu muốn cho các hình gắn vào lưới, ta chọn Snap trong menu Option. 
Nếu muốn tỉ lệ xích của trục x và t bằng nhau để hình tròn nhìn không 
giống hình ellip ta chọn Axes Equal trong menu Option. 
) Để vẽ miền Ω  ta dùng menu Draw hay các  icon  trên  thanh công cụ 
ngay phía dưới các menu. 
)   Để đặt điều kiện biên ta dùng menu Boundary hay icon ∂Ω. Ta bấm 
lên từng đoạn biên để đặt điều kiện cho nó. 
 435
)  Tiếp theo ta tạo lưới bằng cách dùng menu Mesh hay icon ∆. Để tinh 
chỉnh lưới ta bấm vào Refine Mesh hay icon  
) Tiếp  theo  ta mô  tả dạng phương  trình và  các  thông  số  của nó bằng 
cách dùng menu PDE. Muốn  thế,  ta mở menu PDE hay chọn  icon PDE 
và chọn PDE Specification và cho các tham số của phương trình. 
)   Để giải phương trình ta dùng menu Solve hay chọn icon =  . Ta chọn 
menu  con Parameters  để nhập  điều kiện  đầu và khoảng  thời gian  tìm 
nghiệm 
) Nếu muốn vẽ kết quả, ta dùng menu Plot 
3. Một số ví dụ: 
  a. Ví dụ 1: 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.1) 
trong miền 0 ≤ x ≤ 4, 0 ≤ y ≤ 4 với các điều kiện biên: 
  u(0, y) = ey ‐ cosy   u(4, y) = eycos4 ‐ e4cosy       (vd1.2) 
u(x, 0) = cosx ‐ ex    u(x, 4) = e4cosx ‐ excos4       (vd1.3) 
Để giải phương trình ta thực hiện các bước sau: 
¾ Mở công cụ PDETOOL. Vào menu Option | Axes Limit để hiệu chỉnh 
lại phạm vi giá  trị  của x và y  là  [0 5]  rồi  chọn Apply và Close. Chọn 
Option | Axes Equal 
¾ Bấm vào icon    để vẽ hình vuông. Khi vẽ xong, nếu chưa đúng kích 
thước ta bấm đúp vào đối tượng bây giờ có tên là R1 để hiệu chỉnh lại 
thành Left: 0, Bottom: 0, Height: 4, Width: 4. 
¾ Bấm vào icon ∂Ω thì đường biên của đối tượng có màu đỏ. Trên mỗi 
đoạn biên ta cho điều kiện biên theo (vd1.2) và (vd1.3). Để ghi điều kiện 
biên cho đoạn nào ta bấm đúp chuột lên đoạn đó. Điều kiện biên đã cho 
là điều kiện biên Dirrichlet. Trên biên trái, ta ghi điều kiện biên: 
h = 1, r = exp(y) ‐ cos(y) 
trên biên phải:  
h = 1, r = eycos4 ‐ e4cosy  
trên biên dưới:  
h = 1, r = cosx ‐ ex  
và trên biên trên:  
h = 1, r = e4cosx ‐ excos4 
 436
¾ Bấm đúp chuột vào icon PDE và chọn phương trình dạng elliptic và 
các thông số theo (vd1.1): c = 1, a = 0, f = 0 
¾ Bấm đúp chuột vào icon      để tạo lưới và sau đó tinh chỉnh nó.  
¾ Bấm đúp chuột vào icon = để giải phương trình. 
¾ Vào menu Plot | Parameters để chọn cách vẽ và sau đó vẽ ra kết quả  
  b. Ví dụ 2: Giải phương trình parabolic: 
  −
⎡ ⎤∂ ∂ ∂+ =⎢ ⎥∂ ∂ ∂⎣ ⎦
2 2
4
2 2
u(x,y,t) u(x,y,t) u(x,y,t)10
x y t
        (vd2.1) 
trong miền 0 ≤ x ≤ 4, 0 ≤ y ≤ 4 và 0 ≤ t ≤ 5000  với các điều kiện đầu và điều 
biên: 
  u(x, y, 0) = 0                  (vd2.2a) 
  u(x, y, t) = eycosx ‐ excosy với x = 0, x = 4, y = 0, y = 4      (vd2.2b) 
Để giải phương trình ta theo các bước sau: 
¾ Mở công cụ PDETOOL. Vào menu Option | Axes Limit để hiệu chỉnh 
lại phạm vi giá  trị  của x và y  là  [0 4]  rồi  chọn Apply và Close. Chọn 
Option | Axes Equal 
¾ Bấm vào icon    để vẽ hình vuông. Khi vẽ xong, nếu chưa đúng kích 
thước ta bấm đúp vào đối tượng bây giờ có tên là R1 để hiệu chỉnh lại 
thành Left: 0, Bottom: 0, Height: 4, Width: 4. 
¾ Bấm vào icon ∂Ω thì đường biên của đối tượng có màu đỏ. Trên mỗi 
đoạn biên ta cho điều kiện biên theo (vd2.2b). Để ghi điều kiện biên cho 
đoạn nào ta bấm đúp chuột lên đoạn đó. Điều kiện biên đã cho là điều 
kiện biên Dirrichlet. Trên biên trái, ta ghi điều kiện biên: 
h = 1, r = exp(y) ‐ cos(y) 
trên biên phải:  
h = 1, r = eycos4 ‐ e4cosy  
trên biên dưới:  
h = 1, r = cosx ‐ ex  
và trên biên trên:  
h = 1, r = e4cosx ‐ excos4 
¾ Bấm đúp chuột vào  icon PDE và chọn phương trình dạng parabolic 
và  các  thông  số  theo  (vd2.1):  c = 1e‐4, a = 0,  f = 0, d = 1. Trong menu 
Solve | Parameters ta ghi Time: 0:100:5000, u(t0) = 0 (điều kiện đầu). 
¾ Bấm đúp chuột vào icon      để tạo lưới và sau đó tinh chỉnh nó.  
¾ Bấm đúp chuột vào icon = để giải phương trình. 
 437
¾ Vào menu Plot | Parameters để chọn cách vẽ và sau đó vẽ ra kết quả  
c. Ví dụ 3: Giải phương trình hyperbolic: 
⎡ ⎤∂ ∂ ∂+ =⎢ ⎥∂ ∂ ∂⎣ ⎦
2 2 2
2 2 2
u(x,y,t) u(x,y,t) u(x,y,t)1
4 x y t
          (vd3.1) 
trong miền 0 ≤ x ≤ 2, 0 ≤ y ≤ 2 và 0 ≤ t ≤ 2  với các điều kiện biên zero và điều 
kiện đầu: 
u(0, y, t) = 0  u(2, y, t) = 0  u(x, 0, t) = 0  u(0, 2, t) = 0  (vd3.2) 
u(x, y, 0) = 0.1sin(πx)sin(πy/2)  ∂u/∂t(x, y, 0) = 0 với t = 0    (vd3.3) 
Để giải phương trình ta theo các bước sau: 
¾ Mở công cụ PDETOOL. Vào menu Option | Axes Limit để hiệu chỉnh 
lại phạm vi giá  trị  của x và y  là  [0 2]  rồi  chọn Apply và Close. Chọn 
Option | Axes Equal 
¾ Bấm vào icon     để vẽ hình vuông. Khi vẽ xong, nếu chưa đúng kích 
thước ta bấm đúp vào đối tượng bây giờ có tên là R1 để hiệu chỉnh lại 
thành Left: 0, Bottom: 0, Height: 2, Width: 2. 
¾ Bấm vào icon ∂Ω thì đường biên của đối tượng có màu đỏ. Trên mỗi 
đoạn biên ta cho điều kiện biên theo (vd3.2). Để ghi điều kiện biên cho 
đoạn nào ta bấm đúp chuột lên đoạn đó. Điều kiện biên đã cho là điều 
kiện biên Dirrichlet. Trên biên trái, ta ghi điều kiện biên: 
h = 1, r = 0 
trên biên phải:  
h = 1, r = 0  
trên biên dưới:  
h = 1, r = 0  
và trên biên trên:  
h = 1, r = 0 
¾ Bấm đúp chuột vào  icon PDE và chọn phương trình dạng parabolic 
và các thông số theo (vd2.1): c = 1/4, a = 0, f = 0, d = 1. Trong menu Solve 
| Parameters ta ghi Time: 0: 0.1: 2, u(t0) = 0.1*sin(pi*x).*sin(pi*y/2). 
¾ Bấm đúp chuột vào icon      để tạo lưới và sau đó tinh chỉnh nó.  
¾ Bấm đúp chuột vào icon = để giải phương trình. 
¾ Vào menu Plot | Parameters để chọn cách vẽ và sau đó vẽ ra kết quả