Matlab偏微分方程快速上手:使用pde有限元工具箱求解二维偏微分方程

笔记 专栏收录该内容
24 篇文章 2 订阅

注:本人使用MatlabR2020a版本。

1.pdetoolbox的调用

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

2.绘制定解区域(解的定义域)

由图形界面可知,解的定义域是 x , y x,y x,y二维坐标构成的平面空间。我们必须设置自己的定解区域,才能定义自己的方程:
导航栏下方的前5个按钮,分别对应绘制矩形求解区域、绘制按中心生成的矩形求解区域、绘制椭圆形(圆形)定解区域、绘制按中心生成的椭圆形(圆形)定解区域、绘制多边形求解区域。使用时,只需要点击后在绘图区域拖拽(多边形除外,多边形区域是在绘图区域点点以确定顶点),就可以生成定解区域了。
这是前五个按钮
上面这是前五个按钮。
在这里,作者随意绘制了一个椭圆形区域,和一个矩形区域。操作的时候用鼠标拖动操作柄拖拽就可以。(也可以画好几个叠起来)
在这里,作者随意绘制了一个椭圆形区域,和一个矩形区域
真的是“随便”画一个就可以哦,因为在matlab下,求解区域的位置坐标精度达到了 1 0 − 16 10^{-16} 1016左右,手动画几乎不可能画准。所以下一步教大家怎么细致地调节边界的坐标。

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,nu=0.
如果用了Dirichlet边界条件,边界将显示为红色;Neumann、Robin边界条件将显示蓝色,效果如下。

边界条件

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

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

在这里插入图片描述

工具箱提供的方程通式如下:
1.椭圆型Elliptic,通用数学形式为 − ∇ ⋅ ( c ∇ u ) + a u = f ; -\nabla \cdot(c\nabla u)+au=f ; (cu)+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 ; dtu(cu)+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 ; d2t2u(cu)+au=f;
4.特征值方程Eigenmodes,若 λ \lambda λ为特征值,则数学形式为 − ∇ ⋅ ( c ∇ u ) + a u = λ d u . -\nabla \cdot(c\nabla u)+au=\lambda d u . (cu)+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时的结果,如下图所示。
在这里插入图片描述

这就是matlab pdetool工具箱的主要使用方法,本人也是小白一枚,所以欢迎大家批评指正,可以在评论区留下你的想法。

参考:偏微分方程(姜礼尚《数学物理方程讲义》第三版)(更新完毕,附课件)来自于西北大学数统学院的马老师的数理方程视频课,是按数学系的讲法讲的数理方程,里面有那么两三个视频是讲如何用Matlab求解偏微分方程,如果你懂偏微分方程的话进去听一遍就会了。
感觉应该是疫情期间这位老师的网课视频?那么西北大学的学生们也太幸福了,因为马老师讲的真的很好!人也很好玩哈哈,还去b站里别的老师的微分几何课程下面评论,正巧那个被评论的老师的同学就在b站讲拓扑哈哈,那个老师讲的也特别细我还给听完了(浙江理工庄老师)。扯远了!但是还是安利!!!

  • 13
    点赞
  • 7
    评论
  • 64
    收藏
  • 一键三连
    一键三连
  • 扫一扫,分享海报

©️2021 CSDN 皮肤主题: 数字20 设计师:CSDN官方博客 返回首页
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、C币套餐、付费专栏及课程。

余额充值