偏微分方程的数值解(六)： 偏微分方程的 pdetool 解法

浅夏110

1 图形界面解法简介
对于一般的区域，任意边界条件的偏微分方程，我们可以利用 MATLAB 中 pdetool 提供的偏微分方程用户图形界面解法。 图形界面解法步骤大致上为：
0 S: C% ~( k6 o4 P$ i
9 y2 x: |8 N% B9 z（1）定义 PDE 问题，包括二维空间范围，边界条件以及 PDE 系数等。
6 S" }; M2 B' T. P; F/ r8 Z- N5 `4 X
, W) a/ w" C" M3 w% s: t（2）产生离散化之点，并将原 PDE 方程式离散化。

（3）利用有限元素法(finite element method；FEM)求解并显示答案。
( h0 X5 {% V6 h2 `% H
- C" U5 G0 w/ \6 g6 Q0 c在说明此解法工具之前，先介绍此 PDE 图形界面的菜单下方的功能图标(icon)按 钮。


5 h4 F9 L' `# E! ~, T0 _9 T; a
2 图形界面解法的使用步骤
要利用 pdetool 接口求解之前，需先定义 PDE 问题，其包含三大部份：
" d6 F( P! j% C. M
（1）利用绘图(draw)模式，定义 需要求解的问题的空间范围(domain)Ω 。
  j; C/ w+ U/ _  _& x4 j6 ]* k
（2）利用 boundary 模式，指定边界条件。
7 d; Q( E* m# [" M' n' c# X9 I) R
3 m4 @; R" C2 c% }5 E. `（3）利用 PDE 模式，指定 PDE 系数，即输入 c，a，f 和 d 等 PDE 模式中的系数。
' W" _' E; ]. y& x1 z4 J  j
在定义 PDE 问题之后，可依以下两个步骤求解
) Z* ?) M, ]) P/ ~2 E# o
（1）在 mesh 模式下，产生 mesh 点，以便将原问题离散化。

（2）在 solve 模式下，求解。

（3）最后，在 Plot 模式下，显示答案。
4 N7 x# [  Y" z* N+ I$ E
2 `; w5 W* g) ]! n% _

& z* g, l5 I2 D0 A
, M- v- W8 o* _  ]
/ z- U6 s; ^% I0 R& u
, {( v2 f9 |) F$ M+ |, m8 l

9 Q  F! h; M" {- E! Z$ e1 m7 {$ b
  g$ k' `3 _1 V2 s2 y5 @% `  L
注意：
7 P% q% {- I0 G3 b) B  L+ e
1. MATLAB 会以图形的方式展示结果，使用者亦可点选 plot 下之“parameters”功 能，选择适当的方式显示图形及数据。例如用 3D 方式显示求解结果。参数设置见图 10， 显示结果见图 11。

& |9 W  n: U7 ?# I  s* i( N9 ~& `
) _0 G  {  s7 ^0 S" M) y' Z& e' k/ H  X
2. 另外，若使用者欲将结果输出到命令窗口中，以供后续处理，可利用 solve 功能 项下之“export solution”指定变量名称来完成。

3. 如果求抛物型或双曲型方程的数值解，还需要通过“solve”菜单下的 “parameters…”选项设置初值条件。
9 x3 |1 J9 L- e4 Z1 x
# w$ G( r' c- s; R: j* e4. 在上面定义边界条件和初始条件时，可以使用一些内置变量。
5 u2 h- s. i* S7 D4 P& D1 G# O
% E/ s8 p+ h) v8 v（1）在边界条件输入框中，可以使用如下变量： 二维坐标 x 和 y，边界线段长度参数（s s 是以箭头的方向沿边界线段从 0 增加到 1）， 外法向矢量的分量 nx 和 ny（如果需要边界的切线方向，可以通过 tx=-ny 和 ty=nx 表示）， 解 u。
# B! Y* B3 F, Z) `  A$ H
/ K& n0 [3 X8 w8 `（2）在初值条件的输入框中，也可以输入用户定义的 MATLAB 可接受变量（p， e，t，x，y）的函数。
$ u; M  G* |0 n9 }! s+ H) `6 ^
例 11 使用 PDETOOL 重新求例 8 的数值解。
- ?" k6 h5 R# i+ Y! c
; ]( c$ k: U4 x- t& q  c                例8  求解正方形区域{(x, y) | −1 ≤ x, y ≤ 1}上的热传导方程
$ M! {/ ~: ?. d) [


边界条件为Dirichlet条件u = 0。

" l+ C. @" f, z+ h& V, U解 这里是抛物型方程，其中c = 1, a = 0, f = 0, d = 1。

1）定义 PDE 问题，包括二维空间范围，边界条件以及 PDE 系数等。我们这里就 省略了。
$ N  n8 ^8 e' a3 Y* `
, l3 `# ~. |" B, A# c& Z# m2）区域剖分以后，通过“Mesh”菜单下的“Export Mesh…”选项可以把 p，e，t 三个参数分别输出到工作间。
# c1 A6 W) F" d* F2 [5 h+ @
/ |+ D. F  D% ^' D9 @3）然后编写函数 fun1(x,y)如下：

function f=fun1(x,y);
f=zeros(length(x),1);
ix=find(x.^2+y.^2<0.16);
0 z- U+ n* _8 |0 E' U( M9 z5 ~f(ix)=1;

6 K3 _$ S, k9 Z4 j其中的变量 x，y 是 MATLAB 可接受的内置变量。 设置“solve”菜单下的“parameters…”选项如下：

: G. X( C+ l7 l& |+ \# ?5 c6 Z! `时间框中输入：linspace(0,0.1,20)；

初值框中输入：fun1。
/ C; X* {+ u9 X* g
4）设置“plot"菜单下的“parameters…”选项如下：选择 Height(3-D plot) 和 Animation 两项。
* W' k- V$ Y4 C  D3 m& B
5）用鼠标点一下工具栏上的“＝”按钮，就可以画出数值解的 3-D 图形。
2 R" z$ X1 C. O  ^. r) f3 |
* z+ e- y( m( l% r5 |0 J1 j' ~
/ ^& Z$ w. L4 Z. g" ]————————————————
! j  d( c0 K! `9 N9 j8 A* |版权声明：本文为CSDN博主「wamg潇潇」的原创文章，遵循CC 4.0 BY-SA版权协议，转载请附上原文出处链接及本声明。
原文链接：https://blog.csdn.net/qq_29831163/article/details/89712663
0 ]2 S/ |, ^5 v; e! L/ n9 `
0 Q- k9 Y+ D. E/ k; u
/ g6 N% ~/ H) s4 b5 l* n6 y( {6 \