偏微分方程数值解.docx
- 文档编号:9113833
- 上传时间:2023-02-03
- 格式:DOCX
- 页数:41
- 大小:633.50KB
偏微分方程数值解.docx
《偏微分方程数值解.docx》由会员分享,可在线阅读,更多相关《偏微分方程数值解.docx(41页珍藏版)》请在冰豆网上搜索。
偏微分方程数值解
第一章概述
1.1偏微分方程工具箱的功能
偏微分方程工具箱(PDEToolbox)提供了研究和求解空间二维偏微分方程问题的一个强大而又灵活实用的环境。
PDEToolbox的功能包括:
(1)设置PDE(偏微分方程)定解问题,即设置二维定解区域、边界条件以及方程的形式和系数;
(2)用有限元法(FEM)求解PDE数值解;
(3)解的可视化。
无论是高级研究人员还是初学者,在使用PDEToo1box时都会感到非常方便。
只要PDE定解问题的提法正确,那么,启动MATLAB后,在MATLAB工作空间的命令行中键人pdetool,系统立即产生偏微分方程工具箱(PDEToolbox)的图形用户界面(GraphicalUserInterface,简记为GUI),即PDE解的图形环境,这时就可以在它上面画出定解区域、设置方程和边界条件、作网格剖分、求解、作图等工作,详见1.4节中的例子。
我们将在第二章详细介绍GUI的使用,在第二章给出大量典型例子和应用实例。
除了用GUI求解PDE外,也可以用M文件的编程计算更为复杂的问题,详见第三章和第四章
的内容。
1.2PDEToolbox求解的问题及其背景
1.2.1方程类型
PDEToolbox求解的基本方程有椭圆型方程、抛物型方程、双曲型方程、特征值方程、椭圆型方程组以及非线性椭圆型方程。
椭圆型方程:
,
椭圆型方程:
其中
是平面有界区域,c,a,f以及未知数u是定义在
上的实(或复)函数。
抛物型方程:
双曲型方程:
.
特征值方程:
其中d是定义在
上的复函数,
是待求特征值。
在抛物型方程和双曲型方程中,系数c,a,f和d可以依赖于时间t。
可以求解非线性椭圆型方程:
其中c,a,f可以是未知函数u的函数。
还可以求解如下PDE方程组;
利用命令行可以求解高阶方程组。
对于椭圆型方程,可以用自适应网格算法,还能与非线性解结合起来使用。
另外,对于Poission方程还有一个矩形网格的快速求解器。
1.2.2边界条件
(1)Dirichlet条件:
(2)Neumann条件:
其中
是
的边界
上的单位外法向量,
和
是定义在
上的函数。
对于特征值问题仅限于齐次条件:
和
。
对于非线性情形.系数
和
可以依赖于u;对于抛物型方程和双曲型方程,系数可以依赖于时间t。
对于方程组情形,边界条件为
(1)Dirichlet条件:
(2)Neumann条件:
(3)混合边界条件为:
其中
的计算要使得Dirichlet条件满足。
在有限元法中,Dirichlet条件也称为本质边界条件,Neumann条件称为自然边界条件。
1.3如何使用FDEToolbox
1.3.1定解问题的设置
员简单的办法是在PDETool上直接使用图形用户界面(GUl)。
设置定解问题包括三个步骤:
(1)Draw模式:
使用CSG(几何结构实体模型)对话框画几何区域,包括矩形、圆、椭圆和多边形,也可以将它们组合使用。
(2)Boundary模式:
在各个边界段上给出边界条件,
(3)PDE模式:
确定方程的类型、系数c,a,f和dc。
也能够在不同子区域上设置不同的系数(反映材料的性质)。
1.3.2解PDE问题
用GUI解PDE问题主要经过下面两个过程(模式)
(1)Mesh模式;生成网格.自动控制网格参数。
(2)Solve模式:
对于椭圆型方程还能求非线性和自适应解。
对于抛物型和双曲型力程.设置初始边值条件后能求出给定t时刻的解。
对于特征值问题,能求出给定区间内的特征值;求解后可以加密网格再求解。
1.3.3使用Toolbox求解非标准的问题
对于非标准的问题。
可以用PDEToo1box的函数。
或者用FEM(有限元法)求解更为复杂的问题。
1.3.4计算结果的可视化
从GUI能够使用Plot模式实现可视化。
可以使用Color,Height和Vector等作图。
对于抛物型和双曲型方程,还可以生成解的动画。
这些操作通过命令行都很容易实现。
1.3.5应用领域
在应用界面提供了丁如下应用领域
.结构力学——平面应力问题
.结构力学——平面应变问题
.静电场问题
.静磁场问题
.交流电磁场问题
.直流导体介质问题
.热传导问题
.9‘散问题
这些界面都有对话框,它包括PDE的系数、边界条件、解的性质等。
1.4解偏微分方程的一个例子
解Poisson方程
,边界条件为齐次Dirichlet类型。
第一步:
启动MATLABl,键入pdetool,按回车键确定便可启动GUI,然后在Options菜单下选择Grid命令,打开栅格,栅格的使用,能使用户容易确定所绘图形的大小,如图1—1
1--1
第二步:
分步完成平面几何造型:
R1-C1-E1+R2+C2。
用菜单或快捷工具,分别画矩形R1、矩形R2、椭圆E1、圆C1、圆C2。
画圆时,首先选中椭圆工具,按鼠标右键并拖动即可、或者在按ctrI的同时,拖动鼠标也可绘制圆。
然后在Setformula栏,进行编辑并用算术运算将将图形对象名称连接起来,删除默认的表达式键入R1-C1-E1+R2+C2,按等号健得到所需图形。
若需要,还可进行储存.
形成M文件。
选择Boundary菜单中BoundaryMode命令,进入边界模式。
单击Boundary菜单中RemoveA11SubdomainBorders选项,去除子域边界。
如果想将几何信息和边界信息进行存储,应选择Boundary菜单中的ExPortDecomposedGeometry.BoundaryCond’s…命令,将它们分别储存于g,b变量中,通过MATLAB形成M文件。
第三步:
选取边界.单击Boundary菜单中SpecifyBounddyConditions…选项,打开Boundaryconditlons对话框,输入边界条件,如图1—4。
本例取缺省条件。
即将全部边界设为齐次Dirichlet条件,边界颜色显示为红色。
第四步:
选择PDE菜单中PDEMode命令,进入PDE模式。
单击PDE菜单中PDESpecification…选项,打开PDE对话框,设置方程类型。
本例取缺省设置,类型为椭圆型,参数c,a,f分别为1,0,10。
第五步:
选择Mesh菜单中InitializeMesh命令,进行网格剖分。
第六步:
选择Mesh菜单中RefineMesh命令,对网格加密。
第七步:
选择Solve菜单中So1vePDE命令,解偏微分方程并显示图形解。
第八步:
单击Plot菜单中Parameters…选项,打开Plotselection对话框,选中Color,Height(3—DPlot)和Showmesh三项。
然后单击Plot按钮,显示三维图形解。
第九步:
如果要画等值线图和矢量场图,单击Plot菜单中Parameters…选项,打开PlotSelection对话框.选中Contour和Arrows两项。
然后单击P1ot按钮,可显示解的等值线图和矢量场图。
第二章PDE图形用户界面
2.1PDEToolbox菜单
File菜单(如图1-1)
图1-1
New新建一个几何结构实体模型(ConstructiveSolidGeomery,简记为CSG),默认文件名为“Untitled”。
Open…从硬盘装载M文件
Save将在GUI内完成的成果储存到一个M文件中。
SaveAs…将在GUI内完成的成果储存到另外一个M文件中。
Print…将PDE工具箱完成的图形送到打印机内进行硬拷贝。
Exit退出PDE工具图形用户界面。
2Edit菜单(如图1-2)
图1-2
Undo在绘制多边形时退回到上一步操作。
Cut将已选实体剪切到剪贴板上。
Copy将已选实体拷贝到剪贴板上。
Paste…将剪贴板上的实体粘贴到当前几何结构实体模型中。
Clear删除已选的实体。
SelectAll选择当前几何结构实体造型CSG中的所有实体及其边界和字域。
3Options菜单(如图1-3)
图1-3
Grid绘图时打开或关闭栅格。
GridSpacing…调整栅格的大小。
Snap打开或关闭捕捉栅格功能。
AxesLimits…设置绘图轴的坐标范围。
AxesEqual打开或关闭绘图方轴。
TurnoffToolbarHelp关闭工具栏按钮的帮助信息。
Zoom打开或关闭图形缩放功能。
Application选择应用的模式。
Refresh重新显示PDE工具箱中的图形实体。
4Draw菜单(如图1-4)
图1-4
DrawMode进入绘图模式。
Rectangle/square以角点方式画矩形/方行(Ctrl+鼠标)。
Rectangle/square(centered)
以中心方式画矩形/方行(Ctrl+鼠标)。
Ellipse/circle以矩阵角点方式画椭圆/圆(Ctrl+鼠标)。
Ellipse/circle(centered)以中心方式画椭圆/圆(Ctrl+鼠标)。
Polygon画多边形,单击鼠标右键可封闭多边形。
Rotate…旋转已选的图形。
ExportGeometryDescription,SetFormula,Labels…
将几何描述矩阵gd、公式设置字符sf和标识空间矩阵ns输出到主工作空间去。
单击Draw菜单中Rotate…选项,可打开Rotate比对活框,通过输入旋转的角度,可使选择的物体按输入的角度逆时针旋转。
旋转中心的选择如果缺省,则为图形的质心,也可以输入旋转中心坐标。
5Boundary菜单(如图1-5)
图1-5
BoundaryMode进入边界模式。
SpecifyBoundaryConditions…对于已选的边界输入条件,如果没有选择边界,则边界条件适用于所有的边界。
ShowEdgeLabels显示边界区域标识开关,其数据是分解几何矩阵的列数。
ShowSubdomainLabels显示子区域标识开关,其数据是分解几何矩阵中的子域数值。
RemoveSubdomainBorder当图形进行布尔运算时,删除已选取的子域边界。
RemoveAllSubdomainBorders当图形进行布尔运算时,删除所有的子域边界。
ExportDecomposedGeometry,BoundaryCond’s…
将分解几何矩阵g、边界条件矩阵b输出到主工作空间。
选择Boundary菜单中SpecifyBoundaryConditions.命令可定义边界条件。
在打开的Boundarycondition对话框,可对已选的边界输入边界条件。
共有如下三种不同的条件类型:
NeMmann条件这里边界条件是由方程系数q和g确定的,在方程组的情况下(换成方程组模式),q是2ⅹ2矩阵,g是2x1矢量。
Dirichlet条件u定义在边界上,边界条件方程是价h*u=r,这里h是可以选样的权因子(通常为1)。
在方程组情况下,h是2x2矩阵,r是2xl矢量,
混合边界条件(仅适合于方程组情形)它是Dirichlet和Neumann的混合边界条件,q是2x2矩阵,g是2x1矢量,h是1x2矢量,r是一个标量。
6PDE菜单(如图1-6)
图1-6
PDEMode进入偏微分方程模式。
ShowSubdomainLabels显示子区域标识开关。
PDESpecification…调整PDE参数和类型。
ExportPDECoefficients…将当前PDE参数c,a,f,d输出到主工作空间,其参数变量为字符类型。
7Mesh菜单(如图1-7)
图1-7
MeshMode输入网格模式。
InitializeMesh建立和显示初始化三角形网格。
RefineMesh加密当前三角型网格。
JiggleMesh优化网格。
UndoMeshChange退回上一次网格操作。
DisplayTriangleQuality用0~1之间数字化的颜色显示三角形网格的质量,大于0.6的网格可接受的。
ShowNodeLabels显示网格节点标识开关,节点标识数据是点矩阵p的列。
ShowTriangleLabels显示三角形网格标识开关,三角形网格标识数据是三角形矩阵t的列。
Parameters…修改网格生成参数。
ExportMesh输出节点矩阵p、边界矩阵e和三角形矩阵t到主工作空间。
8Solve菜单(如图1-8)
图1-8
SolvePDE对当前的几何结构实体CSG、三角形网格和图形解偏微分方程。
Parameters…调整PDE的参数。
ExportSolution…输出PDE的解矢量u。
如果可行,将计算的特征值1输出到主工作空间。
1.1.9Plot菜单(如图1-9)
图1-9
PlotSolution显示图形解。
Parameters…打开绘图方式对话框。
ExportMovie…如果动画被录制了,则动画矩阵M将输出到主工作空间。
10Window菜单
从Window菜单项,可选择当前打开的所有的MATLAB图形窗口,被选择的窗台移至前台。
11Help菜单
Help…显示帮助信息
About…显示版本信息
1.2PDE工具栏
主菜单下是工具栏,工具栏中喊有许多工具图标按钮,可提供快速、便捷的操作方式。
从左到右5个按钮为绘图模式按钮,紧接着的6个为边界、网格、解方程和图形显示控制功能按钮,最右边的为图形缩放功能键。
(如图1-10)
图1-10
以角点方式画矩形/方行(Ctrl+鼠标)。
以中心方式画矩形/方行(Ctrl+鼠标)。
以矩形角点长轴方式画椭圆/圆(Ctrl+鼠标)。
以中心方式画椭圆/圆(Ctrl+鼠标)。
画多边形,按右键可封闭多边形。
进入边界模式。
打开PDESpecification(偏微分方程类型)对话框。
初始化三角形网格。
加密三角形网格。
解偏微分方程。
打开PlotSelection对话框,确定后给出解的三维图形。
为显示缩放切换按钮。
第三章典型方程及其应用实例
求解PDE问题主要有两种方法,一种是使用图形用户界面,另一种是采用命令行编程。
前者直观简便,而后者更为灵活。
2.1求解椭圆方程的例子
例:
单位圆上的Poisson方程边值问题:
这一问题的精确解为:
若使用图形用户界面(GraphicalUserInterface,简记为GUI),则首先在MATLAB的工作窗口中键入pdetool,按回车键确定,于是出现PDEToolbox窗口。
如果需要坐标网格,单击Options菜单下的Grid选项即可。
下面分步进行操作。
(i)画区域圆单击工具
,大致在(0,0)位置单击鼠标右键同时拖拉鼠标到适当位置松开,绘制圆。
为了保证所绘制的圆是标准的单位圆,在所绘圆上双击,打开ObjectDialog对话框,精确地输入圆心坐标X-center为0、Y-cebter为0及半径Radius为1,然后单击OK按钮,这样单位远已画好。
(ii)设置边界条件单击工具
,图形边界变红,逐段双击边界,打开BoundaryCondition对话框,输入边界条件。
对于同一类型的边界,可以按Shift键,将多个边界同时选择,统一设置边界条件。
本题选择Dirichlet条件,输入h为1,r为0,然后单击OK按钮。
也可以单击Boundary菜单中SpecifyBoundaryConditions…选项,打开BoundaryCondition对话框输入边界条件,如图2-1。
(iii)设置方程单击PDE菜单中PDESpecification…选项,打开PDESpecification对话框,选项方程类型。
本题单击Elliptic,输入c为1,a为0,f为1,然后单击OK按钮,如图2-2。
图2-1
图2-2
(iv)网格剖分单击工具
,或者单击Mesh菜单中InitializeMesh选项,可进行初始网格剖分,这时在PDEToolbox窗口下方的状态栏内显示初始问网格的节点数和三角形单元数。
本题节点数为144个,三角形单元数为254个。
如果需要网格加密,再单击
,或者单击Mesh菜单中RefineMesh选项,这时节点数变为541个,三角形单元数为1016个,如此还可继续加密。
(v)解方程单击工具
,或者单击Solve菜单中Solve菜单中SolvePDE选项,可显示方程色彩解。
如果单击Plot菜单中Parameters…选项,出现PlotSelection对话框,如图2-3,从中可以选择Color,Contour,Arrows,Deformedmesh,Height(3-Dpolt),还可以设置等值线的数目等。
本例中选择Color,Contour,Height(3-Dpolt)和Showmesh四项,然后单击Plot按钮,方程的图形解如图2-4所示。
除了作定解问题解u的图形外,也可以作|gradu|,|cgradu|等图形。
图2-3
图2-4
(vi)与精确解作比较单击Plot菜单中Parameters…选项,打开PlotSelection对话框,在Height(3-Dplot)行Property下拉框中选userentry,且在该行的Userentry输入框中键入u-(1-x.^2-y.^2)/4,单击Plot按钮就可以看到解的绝对误差图形,如图2-5.可见在边界处误差为0。
图2-5
(vii)输出网格节点的编号、单元编号以及节点坐标单击Mesh菜单中ShowNodeLabels选项,再单击工具
或
,即可显示节点编号。
若要输出节点坐标,只需单击Mesh菜单中ExportMesh…选项,这时打开的Export对话框中默认值为p,e,t,这里p,e,t分别表示points(点)、edges(边)、triangles(三角形)数据的变量,单击OK按钮。
然后在MATLAB命令窗口键入p,按回车键确定,即可显示出节点按编号排列的坐标(二维数组);键入e,按回车键,则显示边界线段数据矩阵(7维数组);输入t,按回车键,则显示三角形单元数据矩阵(4维数组)。
(viii)输出解的数值单击Solve菜单中ExportSolution…选项,在打开的Export对话框中输入u,单击OK按钮确定。
再在MATLAB命令窗口中输入u,按回车确定,即显示按节点编号排列的解的数值。
我们也可以用MATLAB程序求解PDE问题,同时显示解的图形:
[p,e,t]=initmesh(‘circleg’,’hmax’,1);
Error=[];err=1;
Whileerr>0.001,
[p,e,t]=refinemesh(‘circleg’,p,e,t);
U=assempde(‘circleb1’,p,e,t,,1,0,1);
Exact=(1-p(1,:
).^2-p(2,:
).^2)’/4;
Err=norm(u-exact,’inf’);
Error=[errorerr];
End
Pdemesh(p,e,t)
Pdesurf(p,t,u)
Pdesurf(p,t,u-exact)
通过命令行键入help+命令函数,如helppdemesh,按回车键,可以调入有关命令函数的定义、参数格式等帮助信息。
2.2求解抛物型方程的例子
例:
考虑一个带有矩形孔的金属板上的热传导问题。
板的左边保持在100
,板的右边热量从板向环境空气定常流动,其他边及内孔边界保持绝缘。
初始
时板的温度为0,于是概括为如下定解问题:
域
的外边界顶点坐标为(-0.5,-0.8),(0.5,-0.8),(0.5,0.8),(-0.5,0.8)。
内边界顶点坐标为(-0.005,-0.4),(0.05,-0.4),(0.05,0.4),(-0.05,0.4)。
使用GUI求解这一问题。
在PDEToolbox窗口的工具栏中选择GenericScalar模式。
(i)区域设置单击
工具,在窗口拖拉出一个矩形,双击矩形区域,在ObjectDialog对话框中输入Left为-0.5,Bottom为-0.8,Width为1,Height为1.6,单击OK按钮,显示矩形区域R1。
用同样方法作内孔R2,只要设置Left=-0.05,Bottom=-0.4,Width=0.1,Height=0.8即可。
然后在Setformula栏中键入R1-R2。
(ii)设置边界条件单击
,使边界变红色,然后分别双击每段边界,打开BoundaryCvondition对话框,设置边界条件。
在左边界条件。
在左边界上,选择Dirichlet条件,输入h为1,r为100;右边界上,选择Neumann条件,输入g为-1,q为0;其他边界上,选择Neumann条件,输入g为0,q为0。
(iii)设置方程类型单击
,打开PDESpecification对话框,设置方程类型为Parabolic(抛物型),d=1,c=1,a=0,f=0,单击OK按钮。
(iv)网格剖分单击
,或者加密网格,单击
。
(v)初值和误差的设置单击Solve菜单中Parameters…选项,打开SolveParameters对话框,输入Time为0:
5,u(t0)为0,Relativetolerance为0.01,Absolutetolerance为0.001,然后单击OK按钮。
(vi)数值解的输出单击Solve菜单中ExportSolution…选项,在Export对话框中输入u,单击OK按钮。
再在MATLAB命令窗口中键入u,按回车键,这时显示出按节点编号的数值解。
(vii)解的图形单击Plot菜单中Parameters…选项,打开PlotSelection对话框,选Color,Contour,Arrows,单击Plot按钮,窗口显示出Time=5时解的彩色图形,如图2-6。
图2-6
MATLAB程序:
[p,e,t]=initmesh(‘crackg’);
u=parabolic(0,0:
0.5:
5,‘crackb’,p,e,t,1,0,0,1);
pdeplot(p,e,t,‘xydata’,u(:
11),‘mesh’,‘off’,‘colormap’,‘hot’)
2.3求解双曲型方程的例子
例:
考虑如下二维波动方程的定解问题:
用GUI求解。
类似前面的例子,首先作正方形区域:
设置断点坐标为(-1,1),(-1,-1),(1,-1),(1,1)。
在ObjectDialog对话框中输入Left为-1,Bottom为-1,Width为2,Height2,单击OK按钮。
设置边界条件:
左、右边界用Dirichlet条件,输入h为1,r为0;上、下边界用Neumann条件,
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 微分方程 数值