数学建模社区-数学中国
标题:
偏微分方程的数值解(六)： 偏微分方程的 pdetool 解法
[打印本页]
作者:
浅夏110
时间:
2020-6-10 10:34
标题:
偏微分方程的数值解(六)： 偏微分方程的 pdetool 解法
1 图形界面解法简介
) F( H( X/ N+ P1 ^ V
对于一般的区域，任意边界条件的偏微分方程，我们可以利用 MATLAB 中 pdetool 提供的偏微分方程用户图形界面解法。 图形界面解法步骤大致上为：
6 m# Y; {& k% F! t" c. d
% Q# D% l& V5 K: k
（1）定义 PDE 问题，包括二维空间范围，边界条件以及 PDE 系数等。
( _1 Z4 L5 B2 X* m" z
- j$ I+ [- X0 l! B
（2）产生离散化之点，并将原 PDE 方程式离散化。
/ O7 E( h. V* \9 O2 t( g8 ~
0 n5 t0 N: N; y- }! r) P) H% Y
（3）利用有限元素法(finite element method；FEM)求解并显示答案。
6 V3 v2 n7 k; _: H/ s3 h: P
" s2 J. e. M. a6 y
在说明此解法工具之前，先介绍此 PDE 图形界面的菜单下方的功能图标(icon)按 钮。
, p& k& O2 s3 L- k9 ?9 G
( m. o4 ?, [$ z8 \
6 C/ ], p( {* c+ u8 {- [, n+ b
* E3 I) i' F {. c( o
2 图形界面解法的使用步骤
6 t8 V% J/ l, v7 C
要利用 pdetool 接口求解之前，需先定义 PDE 问题，其包含三大部份：
" ~% k C( ~; d# Z2 Y6 W# f R
3 W# }4 n- u, K% R
（1）利用绘图(draw)模式，定义 需要求解的问题的空间范围(domain)Ω 。
* B; t" W3 D" ~- T
9 h0 M! ?+ K+ D+ T. ^8 x- G
（2）利用 boundary 模式，指定边界条件。
! A" S2 C# T& x2 G
* f) u* s% ^7 q8 U t5 c" l+ j; l
（3）利用 PDE 模式，指定 PDE 系数，即输入 c，a，f 和 d 等 PDE 模式中的系数。
' |: i+ u1 V. t
9 g) w% {3 P5 B
在定义 PDE 问题之后，可依以下两个步骤求解
7 n/ z0 S; e' |+ T. i
9 H* F7 v( \ u0 o# r, p0 X s5 u
（1）在 mesh 模式下，产生 mesh 点，以便将原问题离散化。
2 F9 }& r" U1 M, D" d
5 v3 |* v: \6 o) V8 s' G6 t. A
（2）在 solve 模式下，求解。
9 q/ K' B1 j3 ]/ U0 Y
5 `5 Q; Z$ U% Z
（3）最后，在 Plot 模式下，显示答案。
" i( X* {4 Z1 n6 }5 R* b
. O( c$ U; v! m1 W+ f
; a0 S+ J* d0 U% R+ O: f& v# @
) R" N8 [: J, ?
}( d0 l4 G& C6 n* C
3 N1 L, k$ U1 S) P! k
0 q4 P. x8 _8 j- v5 \
/ W9 u/ F: h9 r8 c
3 Z, X8 L% R' u2 c. F# V+ o
: f3 j4 u1 f' z# ^; \9 |7 W
6 {& ]7 T* Z4 w! `
注意：
0 ^3 {( k' J, ?# P
% ?1 h5 N* c- z5 q/ s$ y4 e
1. MATLAB 会以图形的方式展示结果，使用者亦可点选 plot 下之“parameters”功 能，选择适当的方式显示图形及数据。例如用 3D 方式显示求解结果。参数设置见图 10， 显示结果见图 11。
|- Y+ _7 l0 {, m
l' i! U5 t, R, h. ?; @
% h% P1 Z) v- C" m/ I" V: m1 ]
) L: D2 n+ H$ v4 R, G( g+ T& _) y
2. 另外，若使用者欲将结果输出到命令窗口中，以供后续处理，可利用 solve 功能 项下之“export solution”指定变量名称来完成。
: N$ Z! R \8 e/ a- i
8 @% d! E1 m J
3. 如果求抛物型或双曲型方程的数值解，还需要通过“solve”菜单下的 “parameters…”选项设置初值条件。
+ t5 C4 L8 w3 T/ i7 K
" L( P" W6 c: y, b z, j
4. 在上面定义边界条件和初始条件时，可以使用一些内置变量。
8 M4 H3 I" c" C# M0 L% O
- c5 h' e8 h/ [' u ^, Q2 y9 T' ^
（1）在边界条件输入框中，可以使用如下变量： 二维坐标 x 和 y，边界线段长度参数（s s 是以箭头的方向沿边界线段从 0 增加到 1）， 外法向矢量的分量 nx 和 ny（如果需要边界的切线方向，可以通过 tx=-ny 和 ty=nx 表示）， 解 u。
5 ?4 `6 t2 a! `3 e
2 j8 ?3 U6 z4 a
（2）在初值条件的输入框中，也可以输入用户定义的 MATLAB 可接受变量（p， e，t，x，y）的函数。
W$ W, Y* R, u4 I Z
" e( m7 b, l: t W) B
例 11 使用 PDETOOL 重新求例 8 的数值解。
; ~ k; g3 S5 _( X# a: H: @% A
8 m* G4 P4 ]* x9 k: Q8 H
例8 求解正方形区域{(x, y) | −1 ≤ x, y ≤ 1}上的热传导方程
: U& O H9 J# X( R* r' g7 a
+ i3 U$ c! O) q" B' s6 J( g
9 \4 W" ]: m: v1 F, g" ~
' K/ o& f& o; U2 _8 B3 @$ u6 Z* a
边界条件为Dirichlet条件u = 0。
( c" ^" b: @, o. i8 m
" [' B g/ [2 H+ V% O: a3 s, e
解 这里是抛物型方程，其中c = 1, a = 0, f = 0, d = 1。
" N, V0 h2 s! o3 s8 v! g
: m* r t3 L1 z# a
1）定义 PDE 问题，包括二维空间范围，边界条件以及 PDE 系数等。我们这里就 省略了。
7 D+ L+ j8 a% G
6 d9 R1 x* D; a
2）区域剖分以后，通过“Mesh”菜单下的“Export Mesh…”选项可以把 p，e，t 三个参数分别输出到工作间。
* G8 o, G( F. e i( [
8 @- Q' q6 _, ~/ z: ], ^% M
3）然后编写函数 fun1(x,y)如下：
2 I% q9 J; s- g* G' t, j
/ K3 @' |- S4 b* h! @. O. w- g- p
function f=fun1(x,y);
7 }8 f& S8 F" q8 H$ Z5 K7 O. b
f=zeros(length(x),1);
7 ~5 r ]1 O0 E; N. K# f: O
ix=find(x.^2+y.^2<0.16);
, s) H: \( u( ^+ w g
f(ix)=1;
/ K9 s0 q4 J. X( U4 u
" q0 {" M# Y' N4 {7 I" L( R
其中的变量 x，y 是 MATLAB 可接受的内置变量。 设置“solve”菜单下的“parameters…”选项如下：
_. e2 e4 E5 K: G. ^" G1 L
: K# N* Y+ h3 D/ p$ X
时间框中输入：linspace(0,0.1,20)；
4 o0 F( y- w# C. }/ x) Z) d
+ P. i8 l6 h! ~& o5 f$ o" P- V
初值框中输入：fun1。
! J+ u7 e" r/ E4 ~ y1 v0 O( B3 _4 w8 }2 K
" l0 @" D% ^ {$ w4 ~
4）设置“plot"菜单下的“parameters…”选项如下：选择 Height(3-D plot) 和 Animation 两项。
1 q% @! g* H& S7 S1 \) ?" T* E( M2 r* x
4 g/ }0 i7 k% {1 Z$ ?7 J
5）用鼠标点一下工具栏上的“＝”按钮，就可以画出数值解的 3-D 图形。
1 x& o9 P# u9 Q* E$ A( J( L# f& k- |
7 L' W! e* _4 r, M3 H& x: N1 G
# h1 _6 F( E. H8 z
————————————————
* T3 n2 i# U0 }4 v" ^; o1 v+ X
版权声明：本文为CSDN博主「wamg潇潇」的原创文章，遵循CC 4.0 BY-SA版权协议，转载请附上原文出处链接及本声明。
& S2 t! n1 k8 Q/ `# c: M* {
原文链接：https://blog.csdn.net/qq_29831163/article/details/89712663
2 o) ^( M ^3 Y( K8 N$ U
5 ^$ I+ B2 d. _3 k
. H a M+ b, f
欢迎光临 数学建模社区-数学中国 (http://www.madio.net/)
Powered by Discuz! X2.5
6 C/ ], p( {* c+ u8 {- [, n+ b
/ W9 u/ F: h9 r8 c
% h% P1 Z) v- C" m/ I" V: m1 ]
9 \4 W" ]: m: v1 F, g" ~