电子课件 [数学物理方法与仿真(第3版)][杨华军][电子教案(PPT版本)]chapter22

凡是 “%”后的语句为解释语句,MATLAB不执行）
（1）题目定义
g='squareg';
% 定义单位方形区域
b='squareb3';
% 左右零边界条件，顶底零导数边界条件
c=1;a=0;f=0;d=1;
（2）初始的粗糙网格化
[p,e,t]=initmesh('squareg');
（3）初始条件
22.1 用偏微分方程工具箱求解微分方程
直接使用图形用户界面（Graphical User Interface，简记作GUI）求解．
例 22.1.1 解热传导方程 ut u f
边界条件是齐次类型，定解区域自定。
计算机仿真
【解】 第一步：启动MATLAB，键入命令pdetool并回
车，就进入GUI．在Options菜单下选择Grid命 令，打开栅格．栅格使用户容易确定所绘图形的 大小． 第二步：选定定解区域 本题为自定区域 ：自拟定解区域如图22.1所 示：E1-E2+R1-E3．具体用快捷工具分别画椭 圆E1、圆E2、矩形R1、圆E3．然后在Set formula栏中进行编辑并用算术运算符将图形对 象名称连接起来. (或删去默认的表达式，直接键入E1-E2+R1-E3)
2.动画图形显示 为了将所得的解形象地表示出来，还要通过一些动画图形命 令．为了加速绘图，首先把三角形网格转化成矩形网格．调用形 式如下： (1)uxy=tri2grid（p，t，u1，x，y） p、t是描述三角形网格的矩阵，x、y是求解区域中矩形网格的坐 标点（矩阵x、y必须都是递增顺序），u1是各时刻三角形网格中 的解．输出矩阵uxy是用线性插值法在矩形网格点上得出的相应u 值． (2) [uxy,tn,a2,a3]=tri2grid(p,t,u,x,y) uxy、p、t、u、x、y意义同上，tn是格点的指针矩阵，a2、a3是内 插法的系数． (3) uxy=tri2grid(p,t,u,tn,a2,a3) 用此命令之前，应先用一个tri2grid命令得出矩阵tn、a2、a3．用此 方法可以加快速度．
（7）在0～5时间内动画显示
newplot;
%建立新的坐标系
newplot;
M=moviein(n);
umax=max(max(u1));
umin=min(min(u1));
for i=1: n, ...
%注意‘…’符号不可省略
if rem(i,10) == 0, ... %当n是10的整数倍时，在命令窗口打印出相应的数字
%若要显示持续不断的动画，则再加上下面语句：
nfps=5;
movie(M,10,nfps);
动态解图可以直接通过MATLAB仿真程序执行看出，
图22.7 是动态图的某一瞬间的解的分布。
图 22.7 某一瞬时的波动方程的解图
例 22.1.2 求解下列热传导定解问题．
u (2u t x2 u(x, y,t) |x
图 22.4 偏微分方程的图解图
第七步：单击 Plot 菜单中 Parameter 选项，打开 Plot Selection 对话框，选中 Color, Height(3-D plot)和 Show mesh 三项．再单击 Polt 按钮，显示三维图形解，如图 22.5 所示．
图 24.15
图 22.5 偏微分方程的三维图形解
计算机仿真求解的偏微分方程类型分为：
椭圆型方程： (cu) au f
抛物型方程： d u (cu) au f t
双曲型方程： d 2u (cu) au f t 2
特征值问题： (cu) au du
特征值偏微分方程中不含参数 f ．
❖ （2）用M文件编程求解．
❖ 本章首先对可视化方法（GUI）求解作初步介绍， 然后详细介绍用M文件编程解几类基本偏微分方 程，并对典型偏微分的解的静态（或动态）显示曲 线分布进行了讨论。
t 2
求解双曲型方程调用形式如下：
（a、c、d、f 是参数）．
(1) u1=hyperbolic（u0，ut0，tlist，b，p，e，t，c，a，f，d）
其中：参数 c、a、f、d 决定了方程的类型．b 代表求解域的边界条件．b 可以是边界 条件矩阵，也可以是相应的 PDE 边界条件 M 文件名．
(3) [p,e,t]=refinemesh(g,p,e,t)
即是迭代过程，得到更细小的网格，使结果更精确．
其中：u0、ut0（ut 即是 u / t ）是初始条件．
tlist 是 t=0 时刻以后均匀的时间矩阵． hyperbolic 函数返回的是 u 在 tlist 中各个时间点、在区域各三角形网格处的值 u1．u1 的每行是由 p 中相应列上的坐标值所得的函数值，u1 的每一列是矩阵 tlist 中相 应的时间项所对应的函数值．
第八步：若要画等值线图和矢量场图，单击 Plot 菜单中 Parameter 选项，在 Plot selection 对话框中选中 Contour 和 Arrows 两项．然后单击 Plot 按钮，可显示解的等值 线图和矢量场图，如图 22.6 所示。
图 22.6 解的等值线图和矢量场图
22.2 计算机仿真编程求解偏微分方程
故参数 a,b, c, d 分别是 1,0,10,1.
第五步：选择 Mesh 菜单中 Initialize Mesh 命令， 进行网格剖分． 选择 Mesh 菜单中 Refine Mesh 命令，使网格密集化， 如图 22.3．
图 22.3 网格密集化
❖ 第六步： 解偏微分方程并显示图形解
❖ 选择Solve菜单中Solve PDE命令，解偏微分 方程并显示图形解，如图 22.4 所示。
x=p(1,:)';
% 注意坐标向量都是列向量
y=p(2,:)';
u0=atan(sin(pi/2*x));
ut0=2*cos(pi*x).*exp(cos(pi/2*y));
（4）在时间段0～5内的31个点上求解
n=31;
tlist=linspace(0,5,n);
% 在0～5之间产生n个均匀的时间点
注意:如果矩形网格点在三角形网格之外，则结果中将会出现出错信息 ‘NAN’．主要的绘图（包括动画）命令函数有：moviein、movie、pedplot、 pdesurf等，其具体应用见下面的例题．
例24.2.1 用 MATLAB 求解下面波动方程定解问题并动态显示解的分布
2u t 2
(
2u x2
网格坐标描述矩阵 p,e,t 是由网格初始化命令得到的：
(2) [p,e,t]=initmesh(g)
其中：g 代表求解区域几何形状，是相应的 PDE 几何分类函数 M 文件名． initmesh 函数的作用是将求解区域进行三角形网格化，网格大小由区域的几何形状决 定．输出的 p、e、t 都是网格数据．点阵 p 的第 1、2 行分别包含了网格中点的 x、y 坐标．e 是边缘矩阵，其各行的意义与求解步骤无直接关系，从略．t 是三角矩阵， 其中的几行描述了区域的顶点．有时还会用到修整网格（精细化）命令
2u ) y 2
y1 u
0 |x y1
0
u(
x,
y,
0)
1 0
(r 0.4) (r 0.4)
求解域是方形区域，其中空间坐标的个数由具体问题确定．
图 22.8 某一瞬时的热传导方程解分布
22.1.3 椭圆型：稳定场方程的求解
稳定场方程属于椭圆方程，下面我们来求解含有源项的标准稳定场方程， 也即是泊松方程．在 MATLAB 中是指如下形式：
2u y 2
)
0
u
u
u
|x1
u
|x1
0,
yHale Waihona Puke y1 yy1 0u(x,
y, 0)
atan[sin( π 2
x)], ut
( x,
y, 0)
2 cos(πx) exp[cos( π 2
y)]
已知求解域是方形区域，空间坐标的个数由具体问题确定．
【解】采用步骤如下（注：在MATLAB环境下运行时，下面的中文不输 入。
图22.1 所讨论定解问题的区域
第三步：选取边界 首先选择Boundary菜单中Boundary Mode命
令，进入边界模式．然后单击Boundary菜单中 Remove All Subdomain Borders选项，从而去掉子 域边界，如图22.2．单击Boundary菜单中Specify Boundary Conditions选项，打开Boundary Conditions对话框，输入边界条件．本例取默认条 件，即将全部边界设为齐次Dirichlet条件，边界显 示为红色．如果想将几何与边界信息存储，可选择 Boundary菜单中的Export Decomposed Geometry,Boundary Cond’s命令，将它们分别存储 在g、b变量中，并通过MATLAB形成M文件．
(2) 用有限元法（FEM）求解PDE．即网格的生成、方程的离 散以及求出数值解；
(3) 解的可视化． 用PDE Toolbox可以求解的基本方程有：椭圆方程、抛物方程、 双曲方程、特征值方程、椭圆方程组以及非线性椭圆方程．
具体操作上可用两个途径： (1) 直接使用图形用户界面（Graphical User Interface， 简记作 GUI）求解．
❖ 求解偏微分方程除了上节介绍的直接使用偏微分方程工具 箱外，还可以用编写程序（对于MATLAB仿真，可以编写 M文件）的方法求解偏微分方程。
❖ 22.2.1双曲型：波动方程的求解 ❖ 1. 求解双曲型方程
❖ 下面将讨论标准波动方程的求解问题．波动方程属于双曲 型方程，即
2u d (cu) au f