Bài giảng Matlab - Chương 8: Phương trình vi phân đạo hàm riêng

CHƯƠNG 8: PHƯƠNG TRÌNH VI PHÂN  
ĐẠO HÀM RIÊNG 
§1. MỞ ĐẦU 
1. Khái niệm  chung: Partial Differential Equation  (PDE) Toolbox  cung  cấp 
một môi trường mạnh và mềm mại để nghiên và giải các phương trình vi phân 
đạo hàm riêng trong mặt phẳng. Dạng phương trình cơ bản của PDE Toolbox 
là: 
  ‐∇.(c∇u) + au = f trong miền Ω 
Các phương trình được rời rạc hoá bằng phương pháp phần tử hữu hạn(FEM). 
Các đối tượng trong PDE cung cấp công cụ để: 
• xác định bài toán PDE, nghĩa là xác định vùng 2D, các điều kiện biên và 
các hệ số PDE. 
• giải bằng phương pháp số các bài  toán, nghĩa  là  tạo ra  lưới không có 
cấu trúc, rời rạc hoá phương trình và tìm nghiệm xấp xỉ.  
• hiển thị kết quả 
2. Sử dụng GUI: 
  a. GUI: PDE Toolbox có một bộ giao diện đồ hoạ người dùng (graphical 
user interface GUI) bao gồm các khiá cạnh của quá trình tìm nghiệm của PDE. 
Để kích hoạt nó ta nhập lệnh pdetool tại cửa sổ lệnh của MATLAB. 
  Trên cửa sổ GUI có các menu và các  icon. Ta dùng các menu hay  icon 
này để thực hiện các nhiệm vụ nhất định. 
  b. Các menu: Có một số menu sau đây trên cửa sổ GUI: 
  • File: Từ menu này ta có thể Open và Save mô hình dưới dạng M‐file. Ta 
có thể in đồ thị và thoát khỏi GUI.  
  • Edit:Từ menu này ta có thể cắt, dán, xoá, copy các đối tượng và chọn tất 
cả các đối tượng. 
  • Options: Menu này chứa các tuỳ chọn cũng như các thay đổi trên trục x.  
  • Draw: Từ menu này ta có thể chọn các đối tương cơ bản để vẽ. 
  • Boundary: Menu này dùng để nhập các điều kiện biên cho các vùng. 
  • PDE: Menu này cung cấp các hộp thoại để mô tả PDE và xuất các hệ số 
của nó vào vùng làm việc.  
  • Mesh: Ta dùng menu này để tạo ra lưới và thay đổi các tam giác. 
  • Solve: Ta chọn menu này để giải các phương trình PDE. 
  • Plot: Từ menu này ta vẽ nghiệm. 
154
  • Window: Chọn cửa sổ làm việc 
  • Help: Hiển thị cửa sổ trợ giúp. 
  c. Thanh công cụ: 
Thanh công cụ có dạng như hình vẽ: 
d. Các GUI mode: Quá trình giải PDE gồm các bước sau: 
  • xác định vùng 2‐D 
  • xác định điều kiện biên 
  • xác định PDE 
  • tạo lưới tam giác 
  • giải PDE 
  • vẽ nghiệm và các thuộc tính vật lí 
pdetool GUI được thiết kế theo cách tương tự. Ta  làm việc trong 6 kiểu khác 
nhau, mỗi kiểu tương ứng với một bước trong quá trình giải PDE. 
  • Trong Draw mode ta tạo vùng 2‐D. 
  • Trong Boundary mode ta có thể mô tả điều kiện biên 
  • Trong PDE mode ta nhập các hệ số của phương trình vi phân. 
  • Trong Mesh mode ta khởi tạo lưới. 
  • Trong Solve mode ta giải phương trình. 
  • Trong Plot mode ta ve nghiệm và các thuộc tính vật lí khác 
  e. Mô hình CGS và Set Formular: PDE Toolbox dùng mẫu mô hình CSG 
để mô hình hoá.Ta có thể vẽ các đối tượng chồng nhau. Có 4 loại đối tượng: 
  • Circle 
  • Polygon 
  • Rectangle 
  • Ellipse 
Mỗi đối  tượng có một  tên duy nhất  trong GUI. Tên mặc định của đối  tượng 
circle là C, đối tượng đa giác là P, đối tượng hình chữ nhật là R, đối tượng hình 
vuông  SQ và  đối  tượng  ellip  là E. Tên  được hiển  thị  trên  đối  tượng. Trong 
Draw mode ta có thể thay đổi tên và hình dạng đối tượng bằng cách nhấp đúp 
lên nó. Khi đó một hộp thoại được mở ra và ta thay đổi thông số. Các toán tử + 
, =, * được dùng để kết hợp các vùng. Toán tử + được dùng để tạo tổ hợp, toán 
tử ‐ dùng tạo vùng có phần rỗng và toán tử * dùng để tạo ra vùng có giao nhau 
155
menu con Select All trong menu Edit và chọn điều kiện biên Neumann cho tất 
cả các biên. Sau đó sửa lại điều kiện biên cho hai phía của thanh. Phía trái chọn 
điều  kiện  biên  Dirichlet  với  r  =  100.  Phía  bên  phải  chọn  điều  kiện  biên 
Neumann với g =  ‐10. Bước  tiếp  theo  là mở hộp  thoại PDE Specification và 
nhập vào các hệ số của PDE. PDE parabolic tổng quát mà PDE Toolbox xử lí 
có dạng: 
  fau)u.(
t
ud =+∇∇−∂
∂  
với điều kiện đầu u0 = u(t0) và thời gian tính nghiệm mô tả trong mảng tlist. 
Như vậy trong trường hợp này ta có d = 1, c = 1, a = 0 và f = 0. Khởi gán các 
lưới và làm tinh lại. Điều kiện đầu u0 = 0 và khoảng thời gian được nhập vào 
là [0:0.5:5]. Ta nhập chúng vào hộp thoại Solve Parameters từ menu Solve. Bây 
giờ  ta có  thể giải bài  toán. Để  thấy được quá  trình  truyền nhiệt  ta đánh dấu 
vào ô Animation  trong hộp  thoại Plot selection. Nên chọn màu  là colormap 
hot. Chú ý  là nhiệt độ của khối  tăng rất nhanh. Bài  toán này được  lưu  trong 
ct8_9.m. 
  b.   Phân bố nhiệt trong thanh phóng xạ: Bài  toán phân bố nhiệt này  là 
một ví dụ về bài toán 3‐D PDE parabolic được biến đổi thành bài toán 2‐D nhờ 
dùng  toạ độ  trụ. Ta khảo sát một  thanh phóng xạ hình  trụ. Tại cuối bên  trái 
của thanh nhiệt được gia tăng  liên tục. Đầu cuối bên phải có nhiệt độ không 
đổi. Tại biên bên ngoài, nhiệt được trao đổi với mô trường bằng truyền nhiệt. 
Tại một thời điểm,nhiệt độ được tạo ra không đồng đều trong toàn bộ thanh 
do quá trình phóng xạ. Giả sử ban đầu nhiệt độ bằng 0. Điều này đưa tới bài 
toán sau: 
  f)uk.(
t
uC =∇∇−∂
∂ρ  
Trong đó ρ là mật độ, C là nhiệt dung riêng của thanh, k là hệ số dẫn nhiệt và f 
là nguồn nhiệt phóng xạ. Mật độ của kim loại là 7800kg/m3, nhiệt dung riêng 
là 500Ws/kg0C, độ dẫn nhiệt là 40W/m0C. Nguồn nhiệt là 20000W/m3. Nhiệt độ 
ở một đầu  thanh  là 1000C. Nhiệt độ môi  trường bên ngoài  là 1000C và hệ số 
truyền nhiệt là 50W/m20C. Dòng nhiệt ở cuối bên trái là 5000 W/m2. Nhưng đây 
là bài toán hình trụ, như vậy ta cần biến đổi phương trình, dùng các toạ độ trụ 
r , z và θ. Do tính đối xứng, nghiệm không phụ thuộc θ. Như vậy phương trình 
đã biến đổi là: 
  fr
z
ukr
zr
ukr
rt
uCr =⎟⎠
⎞⎜⎝
⎛
∂
∂
∂
∂−⎟⎠
⎞⎜⎝
⎛
∂
∂
∂
∂−∂
∂ρ  
Điều kiện biên là: 
168
  •  = 5000 ở đầu cuối bên trái của thanh(điều kiện biên Neumann). 
Do điều kiện Neumann tổng quát hoá trong PDE Toolbox là  .(c∇u)+qu = g và 
c phụ thuộc vào r trong bài toán này( c= kr),điều kiện biên này được biểu diễn 
bằng biểu thức 
)uk.(n ∇r
nr
n
r.(c∇u) = 5000r. 
  • u = 100 tại đầu cuối bên phải của thanh(điều kiện biên Dirichlet) 
  •  nr.(k∇u)=  50(100‐u)  tại  biên  bên  ngoài(điều  kiện  biên Neumann  tổng 
quát hoá). Trong PDE Toolbox nó được biểu diễn bằng: 
  nr.(c∇u)+ 50r.u = 50r.100. 
  • trục của hình trụ r = 0 không phải là biên trong bài toán gốc nhưng khi 
biến đổi thành 2‐D thì lại là biên.Ta phải cho một điều kiện biên nr.(c∇u) = 0 tại 
đây. Giá trị đầu là u(t0) = 0 
Mô hình thanh là hình chữ nhật dọc theo trục x và trục y hướng r. Ta vẽ 
hình chữ nhật với các góc (‐1.5,0), (1.5,0), (1.5,0.2), (‐1.5,0.2), nghĩa là cần nhập 
các số [‐1.5  0.0  3  0.2] vào Object Dialog của phần tử R1. Nhập điều kiện biên 
Neumann cho đầu cuối bên trái với q = 0 và g = 5000*y. Nhập điều kiện biên 
Dirichlet cho đầu cuối bên phải với h = 1và r = 100. Đối với biên ngoài dùng 
điều kiện biên Neumann với q = 50*y và g = 50*y*100. Trên trục ta dùng điều 
kiện biên Neumann với q = 0 và g = 0. Các hệ số của phương trình c = 40*y, a = 
0, d = 7800*500*y và f =20000*y.  
3. Các ví dụ về bài toán hyperbolic: 
a. Phương trình sóng: Ta khảo sát sóng tạo ra từ dao động của một màng 
hình vuông  có  các góc  (‐1,‐1),(‐1,1),(1,‐1) và  (1,1). Phương  trình dao  động  có 
dạng: 
0u
t
u
2
2
=∆−∂
∂  
Màng được cố định(u = 0) tại cạnh phải và cạnh trái và tự do ⎟⎠
⎞⎜⎝
⎛ =∂
∂ 0
n
u ở 
cạnh  trên và cạnh dưới. Ngoài ra,  ta cần giá  trị đầu u(t0) và  ⎟⎠
⎞⎜⎝
⎛
∂
∂
t
)t(u 0 . Giá  trị 
đầu  phải  khớp  với  điều  kiện  biên.  Nếu  ta  bắt  đầu  tại  t  = 
0,thì ⎟⎠
⎞⎜⎝
⎛ π= x
2
cosarctan)0(u và 
⎟⎠
⎞⎜⎝
⎛ π
π=∂
∂ y2sine)xsin(3
t
)0(u  là các giá trị đầu thoả mãn 
điều kiện biên. 
  Ta dùng PDE Toolbox với mode Generic Scalar. Vẽ hình chữ nhật với 
các góc như trên, nghĩa là ta phải điền vào Object Dialog các số: [ ‐1  ‐1  2  2]. 
Sau  đó  ta  xác  định  điều  kiện  biên  và  khởi  gán  lưới.  Mở  hộp  thoại  PDE 
169