快速上手Matlab偏微分方程：使用pde工具箱

注：本人使用MatlabR2020a版本。

1.pdetoolbox的调用

打开MatlabR2020a，在命令行键入pdetool，进入pdetoolbox。

2.绘制定解区域（解的定义域）

由图形界面可知，解的定义域是 x , y x,y x,y二维坐标构成的平面空间。我们必须设置自己的定解区域，才能定义自己的方程：
导航栏下方的前5个按钮，分别对应绘制矩形求解区域、绘制按中心生成的矩形求解区域、绘制椭圆形(圆形)定解区域、绘制按中心生成的椭圆形(圆形)定解区域、绘制多边形求解区域。使用时，只需要点击后在绘图区域拖拽(多边形除外，多边形区域是在绘图区域点点以确定顶点)，就可以生成定解区域了。


   上面这是前五个按钮。

   在这里，作者随意绘制了一个椭圆形区域，和一个矩形区域。操作的时候用鼠标拖动操作柄拖拽就可以。（也可以画好几个叠起来）


   真的是“随便”画一个就可以哦，因为在matlab下，求解区域的位置坐标精度达到了 1 0 − 16 10^{-16} 10−16左右，手动画几乎不可能画准。所以下一步教大家怎么细致地调节边界的坐标。

3.手动调整定解区域的大小

双击刚刚绘制好的区域，弹出一个对话框，里面是我们的定解区域的边界坐标信息（注意不全是坐标），我们可以在这里手动调整定解区域的位置（以矩形区域为例）：


   这就是刚刚打开时的样子。因为这个区域是作者随便画的，所以坐标信息就像随机数一样。下面我们输入精确的数值：Left: -1, Bottom: -1, Width: 2, Height: 2, Name 就用默认的就好。

参数的意义：Left：左边界的坐标(x左)，Bottom：底边界的坐标(y底)，Width：区域宽度，Height：区域高度。这样就得到了 x ∈ [ − 1 , 1 ] , y ∈ [ − 1 , 1 ] x\in[-1,1],y\in[-1,1] x∈[−1,1],y∈[−1,1]的矩形求解区域。调整好的求解区域显示效果如下。

4.调整绘图窗口的显示区域（调整显示坐标限）

有时候我们会发现我们的定解区域太大了，绘图窗口显示不下；或者定解区域太小了，看上去非常不协调。上一个例子中作者的纵坐标显然非常吻合，但是横坐标多出来了（显示了左右两边的白框），那么我强烈建议大家调整完定解区域的坐标以后，再调整一下绘图窗口的显示区域。

   点击导航栏Options，再点击Axes Limits…(意为调整坐标限)，可以手动设置坐标限。这里作者勾选了Auto，这样matlab将自动帮我们调整坐标限，使得定解区域位于界面中央。Options还有其他操作，大家可以自己尝试一下，这里就不介绍了。

调整后的效果如下。

5.确定边界条件

点击导航栏下方第6个按钮( ∂ Ω \partial \Omega ∂Ω，意为 Ω \Omega Ω的边界条件)，它用来显示边界。下面仍然以矩形区域为例。

此时显示了4个边界。双击任意一个边界，会弹出边界条件对话框，可以在这里随意设置边界条件。

在这里可以设置Dirichlet边界条件、Neumann边界条件和Robin边界条件（分别是第一、第二、第三边界条件），其中Robin边界条件和Neumann边界条件集成到一起了。按提示输入对应的系数就可。

在这里作者使用了如下的边界条件：
x = ± 1 , u = 0 ; y = ± 1 , ∂ u ∂ n = 0. x=\pm1,u=0;y=\pm1,\frac{\partial u}{\partial n }=0. x=±1,u=0;y=±1,∂n∂u​=0.

   如果用了Dirichlet边界条件，边界将显示为红色；Neumann、Robin边界条件将显示蓝色，效果如下。

6.确定偏微分方程的形式

点击第7个图标（显示PDE字样），按提示输入偏微分方程的系数即可。在这里笔者求解波动方程： ∂ 2 u ∂ 2 t = ∇ u . \frac{\partial^2 u}{\partial^2 t}=\nabla u. ∂2t∂2u​=∇u.

工具箱提供的方程通式如下：

   1.椭圆型Elliptic，通用数学形式为 − ∇ ⋅ ( c ∇ u ) + a u = f ; -\nabla \cdot(c\nabla u)+au=f ; −∇⋅(c∇u)+au=f;

   2.抛物型Parabolic，通用数学形式为 d ∂ u ∂ t − ∇ ⋅ ( c ∇ u ) + a u = f ; d\frac{\partial u}{\partial t}-\nabla \cdot(c\nabla u)+au=f ; d∂t∂u​−∇⋅(c∇u)+au=f;

   3.双曲型Hyperbolic，通用数学形式为 d ∂ 2 u ∂ 2 t − ∇ ⋅ ( c ∇ u ) + a u = f ; d\frac{\partial^2 u}{\partial^2 t}-\nabla \cdot(c\nabla u)+au=f ; d∂2t∂2u​−∇⋅(c∇u)+au=f;

   4.特征值方程Eigenmodes，若 λ \lambda λ为特征值，则数学形式为 − ∇ ⋅ ( c ∇ u ) + a u = λ d u . -\nabla \cdot(c\nabla u)+au=\lambda d u . −∇⋅(c∇u)+au=λdu.

也可以自行指定求解的方程类型，比如比较常见的热传导、扩散等方程，可以在下面图示的下拉菜单中选择，但是仍然要按上面讲的方法手动设置系数。

7.三角剖分

由于Matlab pdetoolbox使用有限元方法求解，所以需要三角剖分。点击第8个图标（1个三角形图样）可以初始化剖分，点击第9个图标（4个三角形图样）可以增加剖分密度，这样可以提高计算精度，但是密度过高内存可能会爆掉，使用要谨慎。

8.设置初始条件，准备求解

点击导航栏Solve，再点击Parameters…，进入求解参数设置器。在这里，第一行Time我们可以设置 t t t的求解范围及步长，默认情况下是不显示步长的（默认显示0:10意为从0求解到10，步长为1），我们按照Maltab等差数列的生成方法 a:j:b 就可以设置时间步长j了。

第二行u(t0)、第三行u’(t0)表示 t 0 t_0 t0​时刻的两个初始条件。这里作者使用了如下的初始条件。

   第四行和第五行表示相对容差和绝对容差，笔者查看了Matlab帮助中心，大概了解到这两个参数似乎与浮点数0的截断精度有关，太小的话会延长计算时间，如果你想了解更多，笔者把链接提供上来Absolute tolerance - MATLAB & Simulink - MathWorks 中国，假如我们对计算精度没有要求的话，使用默认值就可以了。这里笔者为了演示使用了0.001和0.0001。如果想跟着一起做，那么笔者把方程的代码也放上来：第一个是atan(cos(pi/2*x))，第二个是3*sin(pi*x).*exp(cos(pi*y))。

9.求解

点击导航栏下方按钮(一个“ = = =”字样的按钮，就是增加三角剖分密度右边那个按钮)，这个按钮表示开始求解。如果求解完成的话会显示这个图。在这里可以点击“放大镜”按钮寻找感兴趣的区域放大来观察细节（放大之后想要缩小就要用上面步骤4的方法重新设置坐标限了，没有找到缩小的快捷键）。

   它直接显示了 t = 10 t=10 t=10时 u u u在 求 解 区 域 Ω 求解区域\Omega 求解区域Ω的图像。这样的输出缺乏直观性，我们点击导航栏下方一个长得像matlab的logo的按钮（就是“ = = =”按钮右边那个），调整绘图格式。

   这个窗口有许多功能，作者就不再一一详述了。大家可以自行调试。比较常用的有“Contour”绘制等高线图，“Arrows”绘制向量场，“Height(3-D plot)”按3D模式输出（这个比较常用），“Animation”按动画形式输出(2D\3D都支持)，此选项勾选后右边的Option选项会变亮，我们可以点进去在里面设置1秒显示的帧数、重复播放的次数。这里作者按照25阶变色、‘jet’ Colormap、向量场为 − ∇ u -\nabla u −∇u的形式静态输出 t = 0.5 t=0.5 t=0.5时的结果，如下图所示。

免责声明：本文系网络转载或改编，未找到原创作者，版权归原作者所有。如涉及版权，请联系删