Matlab偏微分方程快速上手:使用PDE toolbox求解二维偏微分方程
作者:互联网
注:本人使用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时刻的两个初始条件。这里作者使用了如下的初始条件。
第四行和第五行表示相对容差和绝对容差,如果对计算精度没有要求的话可以用默认值。这里笔者为了演示使用了0.001和0.0001。
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时的结果。
参考:偏微分方程(姜礼尚《数学物理方程讲义》第三版)(更新完毕,附课件)来自于西北大学数统学院的马老师的数理方程视频课,是按数学系的讲法讲的数理方程,里面有那么两三个视频是讲如何用Matlab求解偏微分方程。感觉应该是疫情期间这位老师的网课视频?感觉西北大学的学生也太幸运了,因为马老师讲的真的很好!人也很好玩哈哈还去b站里别的老师的微分几何课程下面评论,正巧那个被评论的老师的同学就在b站讲拓扑哈哈那个老师讲的特别细我还给听完了(浙江理工庄老师)。扯远了但是还是安利!!!
这就是matlab pdetool工具箱的主要使用方法,本人也是小白一枚,所以欢迎大家批评指正,可以在评论区留下你的想法。
标签:partial,求解,nabla,边界条件,区域,Matlab,微分方程,PDE,坐标 来源: https://blog.csdn.net/weixin_47006934/article/details/113524513