偏微分方程数值解的Matlab 实现

【实例简介】
工程中有许多问题可以归结为偏微分方程问题，如弹塑性力学中研究对象（结构、边坡等）内部的应力应变问题、地下水渗流问题等。这些由偏微分方程及边界条件、初始条件等组合成的数学模型，只有在十分特殊的条件下才能求得解析解。因此，在很长一段时间内，人们对于这一类问题是无能为力的。随着计算机技术的发展，各种数值方法应运而生，如有限元法、有限差分法、离散元法、拉格朗日元法等等。利用数值法，可以求得这些问题的数值解。它不是问题的精确解，但可以无限接近精确解。Matlab采用有限元法求解偏微分方程的数值解
偏微吩方程数宽解的 Matlab罢砚 foxdog制作 中输入公式:C1-C2。 3.3定义边界条件 在 Boundary菜单中选择 Boundary mode选项或直接在工具条中单击a按钮,则显示几何模型的外 边界和内边界。如图24所示。 回区 -1; 11l=3,re,命,=|÷正m=:2-1 〔!1〔2〕001〔E 图24定义边界条件 在图中单击边界,则边界线变成黑色,表示它被选中。按住 Shift键,可以连续选择多条边界线段。 选择要定义的边界以后,双击边界或在 Specify Boundary conditions…近项,打开 Boundary Condition 对话框,如图25所示。 吧ndar丁 Condition 瓦[ F IrI dNrs I: ' ll 1(u-I cue 11 IuI-H 图25 Boundary Condition对话框 对话框中各选项的意义分别为: Boundary cond. cquation标签显示边界条件方程。 Condition type在下面的两个单选钮中进行选择,确定边界条件类型。 Neumann单选钮选择此项,将边界条件类型确定为 Neuman条仵。 Dirichlet单选钮缺省时选此项。选择此项时,将边界条件类型确定为 Dirichlet边界条件。 ● Coelficient value栏在该栏中的对应文丕框中输入边界条件公式卬的系数值。 本例中选择缺省设置。 在 Boundary菜单中近择 Show Edge labels选项和 Show Subdomain labels选项,可以在图中显示 边缘标签和子域标签。 3.4定义PDE类型和PDE系数 偏微分方覆憝宽解的 Matlab寥砚 foxdog制作 在PDE菜单中选择 PdE mode选项,图形将显示为PE模式。选择 Show subdomain labels选项,将 在图中显示子域标签。 在工具条中单击r按钮或在PDE菜单中单击 PDE Specification…选项,可以打开PDE Specification对话框,如图26所丌。 E日 +改 二1三 图26 PDE Specification对话框 对话框中各选项的意义分别为 Equation标签该标签显示PDE方程表达式。 Type of pD栏在该栏中进行选择,确定PDE问题的类型。 ▲ Elliptic单选钮为缺省选项。选择此项,设置为椭圆型PD问题。 Parabolic单选钮选摻此项,设置为抛物线型PDE问题。 ▲ Hyperbolic单选钮选择此项,设置为双曲线型PDE问题。 Eigenmodes单选钮选择此项,设置为特征值PE问题。 ● Coefficient栏在该栏中设置PDE方程的系数值。 3.5三角形网格剖分 在工具烂中单击△按钮,或在Mesh菜单中选择 Initialize mesh选项,可以进行研究域的三角形 网格初始化。如图27所示。 四17训n小网 上 [吕 -s「∵斗-1.三卩U..甲亠[E了 11 在工具条中单击按钮或在Mesh菜单中选择 Refine mesh选项,可以对初始网格迸行加密。如图 28所示。利用加密后的网格进行计算,可以获得具更高精度的解。 11|」| 偏微吩方程数宽解的 Matlab罢砚 foxdog制作 在Mesh菜单中选择 Jiggle mesh选项,可以对网格进行微调。选择 Display triangle图29微调网 格 Quality,将显示三角形的质量图。如图29所示。图中,将三角形的最好质量定为1,并用红色表 ,最差质量为0,用兰色表示,二者之间的三角形质量根据质量得分的大小用红色和兰色的过度色表 小。 在Mesh菜单中选择 Show node labels近项,将在图中显示网格节点的编号,如图30所示。选择 Show Subdomain labels选项,将显示各子域的编号。 其司 C 性 勺 h 12111 图30显示子域编号 在Mesh菜单中选择 Parameters…选项,打开 Mesh parameters刘话框,如图31所示。在该对话框 中进行设置,可以明确与网格剖分有关的参数。 Hesh P Ah rate M Jiggle mesh ggm∈ sh parameters Jiggle mod Number of jiggle iterat ons Refinement method oK 图31 Mesh Parameters对话框 对话框叩各控件的意义分别为 Maximun edge siκe文本框在该文本框中输入数值,确定最大边缘大小。 Mesh growth rate文本框在该文本框中输入网格生长率。缺刍值为1.3。 Jiggle mesh核选框迒择此项,微调网格。缺省时选择此项。 Jiggle mode下拉式列表栏在其中进行选择,确定微调模式。可选模式有on、 optimize minimum 和 optimize mean等三项。 Number of jiggle iterations文木框在该文本框口输入微调迭代的次数。 Refinement method下拉式列表框在该控件中进行冼择,确定进行网格加密的方法,包括 regular和 longest两个选项。具体内容参见前面函数部分。 偏微吩方程数宽解的 Matlab罢砚 foxdog制作 3.6PDE求解 在工具条中单击=按钮或在 Solve菜单口选择 Solve pDe选项,可以对前面定义的PDE问题进行求 解。本例的解如图32所示。 斗 1·3111 图32问题的缺省图解(色谱图) 在 Solve菜单中单击 Parameters…选项,打开 Solve parameters对话框,如图33所示。在该对话 框中进行设置,可以确定求解方法和参数。 rIlI Ir H■"l hlurIl'Itri ruri t-=ul lI-'iyle Ncnlliee l0ler3T 1E4 fuyi'nLrI rur=·urk'l≌. I'TilinI lill IIuII_.'Iuu laurI FI:h nr r h J 图33 Solve Parameters对话框 该对话框中各选项的意义分别为 daptive mode核选框选择此项,即选择自适应模式,系统将自适应生成网格并进行求解。 ● Maximum number of triangles文本框在其中输入三角形的最大个数。 Maximum number of refinements文本框在其中输入网格加密的最大次数。 Triangle selection method在下面的三个单选钮中进行选择,确定三角形的选择方法。 Worst triangles单选钮为缺省选项。选择此项,根据最坏三角形进行迒择。 ▲ Relate tolerance单选钮选此项,根据相对容限进行选择。 ▲User- defined function单选钮选择此项,在下面的文本框中输入函数,用该函数选择三角 形 Worst triangle fractions文本框在该文本框中翰入最坏三角形分量。 Refinement method下拉式列表框在其中进行选择,确定进行网格加密的方法,有 regular 和 longest两个选项。 User nonlincar solver核选框选择此项,用非线性求解器进行求解。 偏微吩方程数宽解的 Matlab罢砚 foxdog制作 onl inear tolerance文本框在该文本中输入非线性迭代终止的容限。缺省值为lE-4。 Initial solution文本框在该文本框中输入问题的解的初值。 Jacobian下拉式列表框在该控件中进行选择,确定雅可比矩阵的确定方式。有 fixed、 lumped 和fu11三个选项。 Norm文本框在其口输入范数值,缺省值为Inf(∞)。 3.7解的图形表达 PDE工具项提供了解的多种图形表达方式,上面显示的是彩色图,为缺省显示。单击顽按钮或在 Solve菜单中单击 Parameters…选项,可以打开 Plot selection刈话框,如图34所示。 图34P1 ot selection对话框 F|t/∈ P=:s; tional 厂 AlmIr,Ad ITE: leigl H plot I cor:i I5 T 厂 animato 甲Flls-』 o- SucIm.l TIn 1.nn 该对话框中各选项的意义分别为: Plot type控件列该列控件控制图形类型的选择。包括 ▲ Color核选框选择此项,生成并显示解的彩色图。缺省时选择此项。 Contour核选框选择此项,生成并显示解的等值线图。如图35所示。 Δ rrows核选框选择此项,生成并显示解的矢量图。如图36所示 De formed mesh核选框选择此项,生成并显示解的变形Ⅸ格图。如图37所示。 ▲ Height(3- D plot)核选框选择此项,生成并显示解的三维图。如图38所示。 Animation核选框选择此顶,生成解的系列演示图。 Property空件列该列控件为一组下拉式列表框,定义对解的哪一部分进行图形显小。 User entry控件列为一组文本框,在其中输入用户输入 Plot style控件列该控件列控制前面选择图型的不冋风格。 Plot in x- y grid核选框选择此项,在x-y网格中绘图。 how mesh核选框选择此项,在当前图中显示网格。 Color: u D 0 E 0 4 0 0.3 " 1-08-0.3-14-0.2C320.40.603112141161 偏微吩方程数宽解的 Matlab罢砚 foxdog制作 图35色谱图叠加等值线图 .2 4 U.2 06 n8-3-14∩.2C12自406811214161.8 图36色谱图叠加矢量图 Y K l 4 凵2「2:斗 1叫1 图37变形网格图 E 图38三维高程图 偏微吩方程数宽解的 Matlab罢砚 foxdog制作 Contour plot levels文本框在其中输入等值线的水平数。 Colormap下拉式列表框在该控件中选择绘彩色图的颜色。 Plot solution automatically核选框选择此项,系统自动绘解的图形。 4几种常见的偏微分方程数值求解问题 4.1椭圆型问题 4.1.1单位圆盘的泊松方程 泊松方程是最简单的椭圆型PDE问颎。 该问题的公式为 △U=1 边界上U=0。该问题的精确解为: U/(x,y)= y (1)、用图形用户界面计算 在命令窗囗中输入 pdetool命令,选用 Generic scalar模式 1.单击0 ption菜单,选择 add a grid选项,设置“snap-to-grid”特点。单击+按钮面一个圆, 若该图不是标准的单位圆,双击该圆,打开一对话棰,在其中可以指定圆心的糈确位置和半径的 大小\。 2.通过单击2按钮来设置边界模式。分割的几何边界显示出来,并目外边界指定为缺省设置,即 Dirichlet边界条件,u=0。本例中采月缺省设置,若边界条件不同,可以通过双击边界打开 对话框,在其中入对应的边界条件。 3.单击按丑,定义偏微分方稈,该操作打开一对话框,可以在其中定义PDE系数c,a和f。本 例中,它们均为常数:c-1,f-1,a-0。 4.单击△按扭或选择Mesh菜单中的 Initialize mesh选项,初始化显示三角形Ⅸ格。 5.单击公按扭或在Mesh莫单中选择 Refine mesh选项,改进初始网格并显示新网格。 0.8 0.6 0.4 0.2 0 -0.2 -0.4 -0.6 -0.8 -0.5 0.5 图39初始化网格 偏微吩方程数宽解的 Matlab罢砚 foxdog制作 6.单击=按扭进行求解,Mat1ab可以用图形来表示问题的解。单击按钮,打开 Plot selection 对语框。利月该对话框,可以选择不同类型的解图。 11 [IT 图40问题解的色谱图 7.比较数值解与精确解之间的误差。 为 Plot selecting对话框中的 Color选择 Property弹出式菜单中的 User entry。然后在 User entry 编辑区中输入 Matlab表达式u-(1-x.2-y.2)/4。可以获得解的绝对误差的图形表示。 图41解的绝对误差图形 (二)、使用命令行函数 首先必须创建 MATLAB函数,使二维几何模型参数化。 M文件 circles.m返回单位圆盘边界点的坐标。下面是该文件的内容 nbs=4 if nargin==0 x=nbs;%边界线段个数 retur eI 0000%参数初值 1111%参数终值 10 【实例截图】
【核心代码】