Chương 8: Phương trình vi phân đạo hàm riêng

Chia sẻ: Nguyen Huu Hao Hao | Ngày: | Loại File: PDF | Số trang:0

Thêm vào BST

Báo xấu

335
lượt xem 135
download

Download Vui lòng tải xuống để xem tài liệu đầy đủ

Tham khảo tài liệu 'chương 8: phương trình vi phân đạo hàm riêng', tài liệu phổ thông, ôn thi đh-cđ phục vụ nhu cầu học tập, nghiên cứu và làm việc hiệu quả

Chủ đề:

Bình luận(0) Đăng nhập để gửi bình luận!

Lưu

Nội dung Text: 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
f. Tạo các góc tròn: Một ví dụ về cách dùng các công thức là tạo ra một hình chữ nhật có góc tròn. Ta khởi động GUI và chọn Grid và Snap to grid. Sau đó thay đổi Grid Spacing trên trục x là [‐1.5:0.1:1.5] và trên trục y là [‐ 1:0.1:1]. Chọn Rectangle/square từ menu Draw để vẽ hình chữ nhật bắt đầu tại(‐1,0.5) rộng 2 và cao 1. Để tạo góc tròn ta thêm các hình tròn ở các góc. Mỗi hình tròn có bán kính 0.2 và tâm cách góc đoạn 0.2. Tiếp đến ta vẽ 4 hình vuông mỗi hình có cạnh 0.2 ở 4 góc. Bây giờ ta soạn công thức cho vùng: R1‐(SQ1+SQ2+SQ3+SQ4)+C1+C2+C3+C4 và ghi vào phần trên của GUI. 3. Sử dụng pdrtool: Để bắt đầu ta sử dụng graphic user interface (GUI) pdetool để giải từng bước PDE. Ta sẽ giải phương trình Poisson ‐∆u = f. Miền 2‐D mà ta giải bài toán sẽ rất phức tạp. Điều kiện biên dùng kiểu Neumann và Dirichlet. Trước hết ta gọi MATLAB. Để gọi GUI ta nhập lệnh pdetool từ cửa sổ MATLAB. Để dễ vẽ hình ta chọn Grid và Snap từ menu Option. Bước đầu tiên để giải bài toán là vẽ vùng 2‐D cần giải. GUI cung cấp 4 loại đối tượng cơ bản: đa giác, hình chữ nhật, và hình ellip. Các đối tượng được dùng để tạo ra mô hình Constructive Solid Geometry(CSG). Mỗi đối tượng được gán cho một nhãn và dùng các nhãn này ta tạo ra được một hình phức tạp là tổ hợp của các hình đơn giản trên. Để chọn một đối tượng, ta nhấp chuột lên icon của nó hay chọn nó từ menu Draw. Icon hình có dấu + dùng để vẽ hình từ tâm. Nếu không có dấu + hình được vẽ từ góc. Khi muốn vẽ hình ta bấm vào icon và đặt con trỏ chuột lên hình rồi kéo và thả. Muốn tạo hình vuông hay tròn ta dùng chuột phải. Đầu tiên ta vẽ một hình chữ nhật và đặt tên là R1. Muốn di chuyển hình ta bấm chuột vào hình và kéo đến vị trí mong muốn rồi thả. Muốn thay đổi kích thước hình bấm đúp chuột vào nó để kích hoạt hộp thoại rồi nhập kích thước mong muốn. Từ hộp thoại ta có thể thay đổi tên của nhãn. Nếu ta muốn xoá hình, chọn nó rồi nhấn Del hay chọn Clear từ menu Edit. Khi được chọn biên đối tượng có màu đen. Muốn chọn nhiều đối tượng cùng lúc ta nhấn thêm phím Shifft. Khi muốn chọn tất cả ta dùng Select All trong menu Edit(hay dùng Ctrl + A). Tiếp theo ta vẽ một hình tròn bằng nhấp chuột vào icon hình ellip có dấu + rồi kéo và thả(dùng chuột phải). CSG ta tạo là hình chữ nhật R1 và E1 có được bằng cách dùng công thức R1 + C1 ở ô Set formula trên cửa sổ. 156
Cuối cùng ta thêm hai đối tượng là hình chữ nhật R2 và hình tròn E2. Để tạo mô hình ta nhập công thức (R1 + E1 + R2) ‐ E2. Ta có thể lưu CSG như là M‐ file bằng cách dùng Save as từ menu File. Sau khi vẽ xong hình ta nhập điều kiện biên cho các biên ngoài cùng. Để nhập điều kiện biên ta bấm chuột vào icon ∂Ω hay chọn Boundary Mode từ menu Boundary. Ta có thể bỏ các biên của vùng con và xác định điều kiện biên. Các đoạn có cạnh xám là biên của vùng con là các vùng chồng nhau của các đối tượng ban đầu. Từ menu Boundary chọn Remove All Subdomain Borders để bỏ tất cả các biên con. Các biên có màu và mũi tên. Màu phản ánh kiểu điều kiện biên và mũi tên hướng về cuối của đoạn biên. Thông tin về hướng được dùng khi khi điều kiện biên được thông số hoá dọc theo biên. Điều kiện biên có thể là hàm của x, y hay là hằng. Mặc định trên biên có điều kiện biên Dirichlet u = 0. Điều kiện biên Dirichlet được thể hiện bằng màu đỏ. Điều kiện biên có thể là Neumann tổng quát(màu xanh) hay hỗn hợp(màu green). Ta chọn điều kiện biên mong muốn bằng cách nhấp chuột vào đoạn biên hay Shift‐click nếu muốn chọn nhiều đoạn biên hay vào menu Edit và chọn Select All và các đoạn biên được chọn có màu đen. Khi này hộp thoại điều kiện biên hiện ra và ta nhập điều kiện biên mong muốn. Trong ví dụ này ∂u ta chọn điều kiện biên Neumann = −5 , nghĩa là độ dốc của nghiệm theo ∂n phương pháp tuyến của đoạn biên này là ‐5. Bấm đúp chuột vào các đoạn biên và trong ô Boundary Condition chọn Neumann. Sau đó nhập giá trị ‐5 vào ô g rồi chọn OK. Giá trị vẫn để mặc định là 0. Sau khi mô tả điều kiện biên ta cần nhập các hệ số của PDE. Ta kích hoạt hộp thoại PDE Specification bằng cách nhấp đúp chuột vào icon PDE hay chọn PDE Specification trong menu PDE. Ta cũng có thể bấm đúp vào từng vùng con để nhập các hệ số PDE cho từng vùng con đó. Trong hộp thoại này ta còn phải chọn kiểu phương trình(elliptic, parabolic, hyperbolic hay eigenmode) và xác định kiểu ứng dụng theo kiểu PDE. Bài toán ta đang xét là bài toán elliptic nên ta đánh dấu vào ô tương ứng. Ta nhập các giá trị c = 1.9, a = 0.0 và f = 10 cho bài toán này. Tiếp theo ta tạo lưới tam giác bằng cách bấm chuột vào icon hình tam giác hay chọn Mesh | Initialize Mesh. Nếu muốn kết quả chính xác hơn ta tinh chỉnh lưới bằng cách nhấp vào icon có 4 tam giác hay chọn Mesh | Refine Mesh. dụng Mesh | Jiggle Mesh ta có thể tăng chất lượng của lưới. Ta có thể huỷ các thay đổi về lưới bằng cách chọn Mesh | Undo. Để giải phương trình ta bầm vào icon = hay chọn Solve | Solve PDE. Kết 157
quả được vẽ ra. Mặc định, hình vẽ dùng màu và thanh màu để chỉ giá trị. Nếu muốn nghiệm được xuất dưới dạng vec tơ cho vùng làm việc của MATLAB. Có nhiều kiểu vẽ cho phép ta quan sát nghiệm. Muốn vậy ta vào menu Plot và chọn kiểu hiển thị nghiệm thích hợp. 4. Dùng các hàm dòng lệnh: Mặc dù pdetool Gui cung cấp một môi trường làm việc thuận tiện nhưng vẫn có những trường hợp ta phải dùng các hàm dòng lệnh. Đó là: • hình dạng hình học của vùng được khảo sát khác đường thẳng, cung tròn, cung ellip. • điều kiện biên không tiêu chuẩn • các hệ số của PDE và của điều kiện biên phức tạp • có trên hai biến phụ thuộc • nghiệm không bị hạn chế • nghiệm khó biểu diễn. Quá trình xác lập bài toán và giải nó được phản ánh trong thiết kế của GUI. Một số cấu trúc dữ liệu xác định các khía cạnh khác nhau của bài toán và các giai đoạn xử lí khác nhau tạo ra các cấu trúc dữ liệu mới. a. Mô hình CSG: Mô hình CSG được mô tả bằng ma trận Geometry Description, set formular và ma trận Name Space. Các cấu trúc dữ liệu này được mô tả trong hàm decsg. b. Phân nhỏ các hình: Hình được phân nhỏ được mô tả bằng ma trận Decomposed Geometry hay bằng ma trận Geometry M‐file. Khi này hình dạng được mô tả như là một bộ các vùng cạnh nhau cực tiểu bao bởi các đoạn biên. Một ma trận Decomposed Geometry có thể được tạo từ mô hình CSG bằng dùng hàm decsg. Nó cũng có thể xuất từ GUI bằng cách chọn Export Decomposed Geometry trong menu Boundary. Một Geometry M‐file tương đương với một ma trận Decomposed Geometry đã cho có thể được tạo ra bằng cách dùng hàm wgeom và xem bằng hàm pdegplot. Cấu trúc dữ liệu của ma trận Decomposed Geometry và M‐file Geometry được mô tả trong hàm pdegeom c. Điều kiện biên: Điều kiện biên được mô tả bằng ma trận Boundary Condition hay M‐file Boundary. Các điều kiện biên được cho như là hàm trên các biên. Ma trận điều kiện biên có thể xuất từ GUI bằng chọn Export Decomposed Geometry, Boundary Cond’s. . .trong menu Boundary. Một M‐file chứa điều kiện biên tương đương với ma trận điều kiện biên có thể tạo được 158
từ hàm wbound. Cấu trúc dữ liệu của ma trận điều kiện biên mô tả trong hàm assemb và pdebound. d. Hệ số của phương trình: PDE được mô tả bằng ma trận Coefficient hay M‐file Coffficient đối với mỗi hệ số c, a, f và d. Các hệ số là hàm trên các vùng con. Các hệ số có thể xuất từ GUI bằng chọn Export PDE Coefficients. . . từ menu PDE. Cấu trúc dữ liệu của nó được mô tả trong hàm assempde. e. Lưới: Một lưới tam giác được mô tả bằng dữ liệu lưới gồm ma trận Points, ma trận Triangle. Trong lưới, vùng nhỏ nhất được tam giác hoá thành các vùng con, các đoạn biên. Số liệu lưới được tạo từ hình dạnh hình học bằng hàm initmesh và thay đổi bằng hàm refinemesh và jigglemesh. Hàm adaptmesh tạo số liệu lưới. Lưới được vẽ bằng hàm pdemesh. f. Nghiệm: Nghiệm của bài toán PDE được biểu diễn bằng vec tơ. Mỗi nghiệm là giá trị tại một điểm lưới của mỗi biến phụ thuộc. Các véc tơ nghiệm được tạo ra bằng hàm assempde, pdenonlin, adaptmesh, parabolic, hyperbolic, và pdeeig. g. Quá trình và biểu diễn: Cho cặp nghiệm‐lưới ta có các công cụ khác nhau để xem các nghiệm. pdeintrp và pdertni có thể dùng để nội suy giữa các hàm. tri2grid dùng để nội suy hàm từ lưới tam giác tới lưới hình chữ nhật. pdegrad và pdecgrad tính gradient của nghiệm. pdeplot để vẽ nghiệm và pdecont, pdesurf thể hiện nghiệm dạng contour và mặt. §1. MỘT SỐ BÀI TOÁN 1. Các ví dụ về bài toán elliptic: a. Phương trình Poisson trên hình tròn đơn vị: Ví dụ đầu tiên về bài toán elliptic là giải phương trình Poisson xác định trên hìng tròn đơn vị. Bài toán được mô tả bằng phương trình: ‐∆u = 1 trong miền Ω, u = 0 trên ∂Ω. Trong đó Ω là hình tròn đơn vị. Trong trường hợp này nghiệm chính xác là: 1 − x2 − y2 u( x , y) = 4 Như vậy ta có thể đánh giá được sai số theo các phương pháp chia lưới khác nhau. Ta thực hiện giải bài toán theo các bước sau: • nhập lệnh pdetool từ cửa sổ lệnh của MATLAB và con trỏ chuột trở thành dấu +. • mở menu Options, đánh dấu mục Grid và Snap. Vẽ hình tròn đơn vị bằng bấm vào icon hình ellip trên thanh công cụ và kéo rồi thả chuột. Nếu 159
hình tròn chưa thỏa mãn yêu cầu thì bấm đúp vào nó để kích hoạt hộp thoại và nhập lại các thông số mô tả chính xác tâm và bán kính của hình tròn. • nhập kiểu biên bằng cách bấm vào menu Boundary và chọn Boundary Mode hay bấm đúp vào nút ∂Ω. Khi này biên của vùng ∂Ω được vẽ và biên ngoài được gán điều kiện biên mặc định(điều kiện biên Dirichlet u = 0 trên biên). Trong trường hợp này đây là điều kiện biên mong đợi. Nếu điều kiện biên khác đi ta bấm đúp vào biên để hiển thị hộp thoại và sửa lại điều kiện biên cho phù hợp và hiển thị nó. • để xác định phương trình vi phân đạo hàm riêng nhấp chuột vào nút PDE trên thanh công cụ(có thể chọn menu PDE | PDE Speficification). Khi này một hộp thoại được mở ra và ta có thể xác định các hệ số c , a và f. Trong trường hợp này c = 1 , f = 1 và a = 0. • nhấp chuột vào nút Mesh(nút hình tam giác) hay chọn Mesh | Initialize Mesh . Khi này một lưới hình tam giác được khởi gán và hiển thị. • nhấp nút Refine(nút tam giác có nhiều tam giác con) hay chọn Mesh | Refine Mesh. Như vậy lưới mịn hơn được khởi gán và hiển thị. • để giải phương trình nhấp nút = hay chọn menu Solve | Solve PDE (Ctrl‐E). Dùng hộp thoại Plot Selection trong menu Plot | Parameters để chọn các hiển thị nghiệm khác nhau. • để so sánh nghiệm số và nghiệm chính xác, chọn menu Plot | Parameters để hiển thị hộp thoại Plot Selection. Trong mục Property của Color chọn user enrty. Sau đó nhập biểu thức MATLAB u‐(1‐x.^2‐y.^)/4 vào trường soạn thảo user etry và nhấp nút lệnh Plot. Ta nhận được hình vẽ sai số tuyệt đối của nghiệm. Để xuất kết quả vào vùng làm việc của MATLAB ta dùng Mesh | Export Mesh và Solve | Export Solution. Để tinh chỉnh kết quả nhấp nút Refine và = nhiều lần. Kết quả lưu trong ct8_2.m. b. Bài toán phản xạ sóng: Bài toán này dùng để tính phản xạ sóng từ một vật thể bị rọi sóng tới. Với bài toán này ta khảo sát một màng nằm ngang rộng vô hạn chịu các dịch chuyển nhỏ theo chiều thẳng đứng u. Màng cố định ở biên. Ta coi môi trường đồng nhất và tốc độ sóng là hằng c. Khi sóng là điều hoà theo t ta có thể tính trường bằng cách giải một bài toán xác lập duy nhất. Với u( x , y , t ) = u( x , y)e − iωt phương trình sóng có dạng: ∂2u − c∆u = 0 ∂t 2 trở thành: − ω2 u − c∆u = 0 160
hay phương trình Helmholz: − ∆u − k 2 u = 0 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 πf 2 π k= = = c c λ 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: r a = (cos(a), sin(a)) : r r v( x , y , t ) = e i( ka .x − ωt ) = v( x , y)e − iωt r r Trong đó : e ika .x = v( x , y) 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: ∂r r + cξ∇r = 0 ∂t 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á: r ξ∇r = ikr 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
r r r = − v( x , y) = −e ika .x 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à: r = − v( x , y) = −e ikx 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: ⎛ ⎞ − ∇.⎜ 1 ∇u ⎟ = 0 ⎜ ⎟ ⎝ 1+ |∇u| 2 ⎠ trong hình tròn đơn vị Ω = {(x , y )|x 2 + y 2 ≤ 1} 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 1 PDE. Đây là phương trình elliptic với c = , a = 0 và f = 0. Hệ số c được 1+ |∇u|2 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: ⎛ K 1 0 ⋅ ⋅ ⋅ 0 B1T ⎞ ⎜ ⎟ ⎜ 0 K 2 ⋅ ⋅ ⋅ 0 B2 ⎟ T K=⎜ M M O M M ⎟ ⎜ ⎟ ⎜ 0 0 L K n BTn ⎟ ⎜ ⎟ ⎝ B1 B 2 L B n C ⎠ Trong khi vế phải là: 163
⎛ f1 ⎞ ⎜ ⎟ ⎜ f2 ⎟ f = ⎜ M ⎟ ⎜ ⎟ ⎜ fn ⎟ ⎜f ⎟ ⎝ c⎠ 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ó: ⎛ K1 0 0 B1T ⎞⎛ u1 ⎞ ⎛ f1 ⎞ ⎜ ⎟⎜ ⎟ ⎜ ⎟ ⎜ 0 K 2 0 B2 ⎟⎜ u 2 ⎟ ⎜ f2 ⎟ T ⎜ 0 = 0 K 3 BT3 ⎟⎜ u 3 ⎟ ⎜ f3 ⎟ ⎜ ⎟⎜ ⎟ ⎜ ⎟ ⎜B B ⎟⎜ c ⎟ ⎜ f ⎟ ⎝ 1 2 B 3 C ⎠⎝ ⎠ ⎝ c ⎠ Và nghiệm xác định bằng cách loại trừ khối: (C − B1K 1−1B1T − B2 K 2−1BT2 − B3K 3−1BT3 )u c = fc − B1K 1−1f1 − B2 K 2−1f2 − B3K 3−1f3 u1 = K 1−1 (f1 − B1T u c ) L 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: ∂u d − ∆u = 0 ∂t 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) ∂u • = −10 ở bên phải (điều kiện biên Neumann) ∂n ∂u • = 0 trên các biên khác(điều kiện biên Neumann) ∂n 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