数学建模社区-数学中国
标题:
偏微分方程的数值解(六)： 偏微分方程的 pdetool 解法
[打印本页]
作者:
浅夏110
时间:
2020-6-10 10:34
标题:
偏微分方程的数值解(六)： 偏微分方程的 pdetool 解法
1 图形界面解法简介
1 R$ @1 B1 a: j7 K
对于一般的区域，任意边界条件的偏微分方程，我们可以利用 MATLAB 中 pdetool 提供的偏微分方程用户图形界面解法。 图形界面解法步骤大致上为：
6 x6 g* u8 R+ I$ y
; c1 o1 d" c! Y6 x
（1）定义 PDE 问题，包括二维空间范围，边界条件以及 PDE 系数等。
: @9 U( ?6 f4 j9 Z, |
! n+ c: f( \. F. a6 w
（2）产生离散化之点，并将原 PDE 方程式离散化。
: Z$ o& B# H, T1 z" e! C
: Y1 w" e, V" m2 I9 Q; }( P; J
（3）利用有限元素法(finite element method；FEM)求解并显示答案。
( F- r2 p* Q) T" G+ n: M2 I
, v/ E/ d% ~# a3 K. T
在说明此解法工具之前，先介绍此 PDE 图形界面的菜单下方的功能图标(icon)按 钮。
1 k2 v9 N$ ?1 q2 h
7 X6 m: s* U7 T! U: H) g% j/ Q
* T# ]+ |" ~& W- G/ O
& B9 ]5 ]! h/ K: q; T; P0 P$ r
2 图形界面解法的使用步骤
! p5 V9 D, E3 T# w1 d/ C8 p {
要利用 pdetool 接口求解之前，需先定义 PDE 问题，其包含三大部份：
/ q( l( L; F8 w. t
# n, _- ^1 H2 r" q
（1）利用绘图(draw)模式，定义 需要求解的问题的空间范围(domain)Ω 。
- X2 m6 k0 B* Z }7 l# y
" g2 |$ W: N( I, b! `
（2）利用 boundary 模式，指定边界条件。
: c+ f0 m1 u5 g. D8 a/ e' x4 w
7 H S) r( k7 K7 Y) R/ k1 M8 z
（3）利用 PDE 模式，指定 PDE 系数，即输入 c，a，f 和 d 等 PDE 模式中的系数。
" _ ^+ g* {- M; G1 }/ e0 [" o
' ^& a( Q* [3 |- E R9 i
在定义 PDE 问题之后，可依以下两个步骤求解
) _% h3 A& Y( y3 D( q6 w
9 [+ h$ `, l; H) C, g
（1）在 mesh 模式下，产生 mesh 点，以便将原问题离散化。
0 K' Y% ]4 k8 }* q% C3 K [1 W0 K
2 A* M4 C! {) _+ f/ |4 o/ W% i" @
（2）在 solve 模式下，求解。
3 a( e0 V9 d# f* Y2 R. d R0 w: V
5 W# o' W4 t/ z4 T
（3）最后，在 Plot 模式下，显示答案。
" y! C m8 x$ u6 d6 p
6 c2 L% `; H' x* _" o; J" J
; o) t& `, Y! ~6 T$ W
* Q0 j8 P/ A# T
; ^9 c: _0 ]! a
8 ?9 N' C1 e3 y1 |4 `; i
9 P* b) E" O* m5 a. N0 R
5 o& S5 r) y0 q" b
* j+ |: F3 C% W L8 B
0 |$ s) G# \# m- H
0 o1 l0 e6 x8 C- e; K
注意：
% R& b( C S [8 @$ c/ h2 v5 r
: T; l l: e" v
1. MATLAB 会以图形的方式展示结果，使用者亦可点选 plot 下之“parameters”功 能，选择适当的方式显示图形及数据。例如用 3D 方式显示求解结果。参数设置见图 10， 显示结果见图 11。
6 p f, B2 ] r+ Q" E# f
" E7 Y' E+ {7 v: J" P( n
& H8 a* [. ?9 l3 h6 M) B( U* L0 u0 Y: ]
3 d6 F% ?: K) t) @7 ^, h
2. 另外，若使用者欲将结果输出到命令窗口中，以供后续处理，可利用 solve 功能 项下之“export solution”指定变量名称来完成。
. U! v- p4 V8 X Z" A' [) X9 E
6 J! T C6 h4 P& ?) P/ v$ Y
3. 如果求抛物型或双曲型方程的数值解，还需要通过“solve”菜单下的 “parameters…”选项设置初值条件。
) e4 o0 M- Q+ I' T
, q% | p+ h9 a1 i9 r7 ~
4. 在上面定义边界条件和初始条件时，可以使用一些内置变量。
' j2 m, `8 w/ V% z: u
1 w* V2 x% \/ A4 z/ Z% a8 X# A: r
（1）在边界条件输入框中，可以使用如下变量： 二维坐标 x 和 y，边界线段长度参数（s s 是以箭头的方向沿边界线段从 0 增加到 1）， 外法向矢量的分量 nx 和 ny（如果需要边界的切线方向，可以通过 tx=-ny 和 ty=nx 表示）， 解 u。
0 n: [( K5 Y7 \$ M, J' M! X
! Y1 o+ y9 F, K2 W0 S2 V
（2）在初值条件的输入框中，也可以输入用户定义的 MATLAB 可接受变量（p， e，t，x，y）的函数。
% J- V, J& f! ?2 l$ ]0 e6 C
+ o% e3 F4 U) w2 Q+ \
例 11 使用 PDETOOL 重新求例 8 的数值解。
- M+ l C; W& Y+ ^1 n, ~! ^+ f
8 |. I" u5 [/ w* ?
例8 求解正方形区域{(x, y) | −1 ≤ x, y ≤ 1}上的热传导方程
' T) H" Q; q+ E' x3 I
" X* n& n9 {; V
2 T, z- a! a3 k4 e
% a$ A) _" q. M0 H' ]" _
边界条件为Dirichlet条件u = 0。
0 V, _% R! P/ b+ W
: Q+ L$ `: A2 G
解 这里是抛物型方程，其中c = 1, a = 0, f = 0, d = 1。
4 A8 a0 _% u4 [5 E
5 l. R; Z! |' C* B+ T
1）定义 PDE 问题，包括二维空间范围，边界条件以及 PDE 系数等。我们这里就 省略了。
. }& o# |6 T0 n! Q, v5 |0 R
, i K. H2 `+ Z! I) h( Q
2）区域剖分以后，通过“Mesh”菜单下的“Export Mesh…”选项可以把 p，e，t 三个参数分别输出到工作间。
% I3 {9 U7 F# ?, S x: s
& k8 m; @) a0 V& U. {
3）然后编写函数 fun1(x,y)如下：
# X% a6 q0 L- Q7 A
% O9 `8 V) N& D
function f=fun1(x,y);
' }2 N# X. Z9 s4 L1 v E+ A
f=zeros(length(x),1);
: l; Q5 n. I. d }8 ?" Y) O2 t% I
ix=find(x.^2+y.^2<0.16);
* ?8 o8 n( m9 B5 m; t4 W
f(ix)=1;
$ f$ x, z- I6 P- `4 \4 L# H
, D( P* f+ g! f" g5 W
其中的变量 x，y 是 MATLAB 可接受的内置变量。 设置“solve”菜单下的“parameters…”选项如下：
2 o' J0 h! P$ D% Y/ C% g% q- O) i, _
+ \- f3 d4 d# \
时间框中输入：linspace(0,0.1,20)；
9 f7 s7 [, r, e& b
3 J6 v; t/ s& F0 I! [5 V( b
初值框中输入：fun1。
; t+ `$ B7 u& s$ i( z
7 H9 R' c. M/ {+ R. o
4）设置“plot"菜单下的“parameters…”选项如下：选择 Height(3-D plot) 和 Animation 两项。
G4 y9 a0 J3 U0 i1 ]+ V
, T+ Z7 H2 Y6 i* F0 @
5）用鼠标点一下工具栏上的“＝”按钮，就可以画出数值解的 3-D 图形。
: j$ i8 B& I& c' n
2 p8 Z: ~! {' U4 ?& r% V$ L
5 Z8 h3 n3 R& Z4 j1 H0 o: t
————————————————
4 K' s; ]+ k2 `
版权声明：本文为CSDN博主「wamg潇潇」的原创文章，遵循CC 4.0 BY-SA版权协议，转载请附上原文出处链接及本声明。
7 L3 [. d% d% G; T
原文链接：https://blog.csdn.net/qq_29831163/article/details/89712663
0 \5 Y/ J7 q' c6 w, A4 G* ?
/ B. r( t& C- f. |% T1 ], K
# Z1 P; i+ S5 V% ~- T
欢迎光临 数学建模社区-数学中国 (http://www.madio.net/)
Powered by Discuz! X2.5
* T# ]+ |" ~& W- G/ O
& H8 a* [. ?9 l3 h6 M) B( U* L0 u0 Y: ]