D用MATLAB解偏微分方程

2006年12月 阴山学刊 Dec.2006 第20卷 第4期 YINSHAN ACADEMIC JOURNAL Vo1.20 No.4 12 用MATLAB解偏微分方程 田 兵 （包头师范学院 学报编辑部，内蒙古 包头 014030） 摘 要：讨论了MATLAB中偏微分方程工具箱的用法。用这个工具箱解方程的过程是：确定待解的偏微分 分方程；确定边界条件；确定方程所在域Ω的几何形状；划分有限元；解方程 关键词：MATLAB；偏微分方程；程序 中图分类号：O175 文献标识码: A 文章编号： 1004－1869（2006）04－0012－02 解偏微分方程不是一件轻松的事情，但是偏微分方程 在自然科学和工程领域应用很广，因此，研究解偏微分方 程的方法、开发解偏微分方程的工具是数学和计算机领域 中的一项重要工作。MATLAB提供了专门用于解二维偏微 分方程的工具箱，使用这个工具箱，一方面解偏微分方 程，另一方面，可以让我们学习如何把求解数学问题的过 程与方法工程化。应当承认，我们国家在数学软件的开发 方面还比较落后，MATLAB是当今世界上最好的数学软件 之一，通过对这个软件的认识，有助于研发我们自己的数 学软件。 MATLAB的偏微分方程工具箱名字叫pdetool，它采用 有限元法解偏微分方程。用这个工具箱可以解如下方程。 椭圆方程 抛物线方程 双曲线方程 特征值方程 所有的方程都在二维平面Ω域上。方程中，▽是 Laplace算子，u是待解的未知函数，c, a, f是已知的实值标 量函数，d是已知的复值函数，λ是未知的特征值。 在边界 上，方程的边界条件一般可以写成 Dirichlet条件(第一类边界条件)：hu=r Neumann条件（第二类边界条件）： 在两个偏微分方程构成方程组的情况下，边界条件可 以写成 Dirichlet条件： Neumann条件 混合条件 式中g, h, q, r是边界 上的复值函数。 是边界 上向外的单位法线。 1 pdetool的使用 在MATLAB命令窗口中键入pdetool，pdetool窗口打开 进入工作状态。pdetool提供两种解方程的方法，一种是通 过函数，利用函数可以编程也可以用命令行的方式解方 程。另一种是对pdetool窗口进行交互操作。一般来说，用 函数解方程比较繁琐，但是比较灵活；通过窗口交互操作 比较简单。解方程的全部过程以及结果都可以输出保存为 文本文件。限于本文的篇幅，我们主要介绍交互操作解偏 微分方程的方法。 1.1 确定待解的偏微分方程 用函数assempde可以对待解的偏微分方程加以描述。 在交互操作中，为了方便用户，pdetool把常见问题归 结为几个类型，可以在pdetool窗口的工具栏上找到选择类 型的弹出菜单。这些类型如下： 通用问题 通用系统（二维的偏微分方程组） ( )+c u au f−∇ ∇ = fauuctu =+∇∇−∂∂ )()/(d fauuctu =+∇∇−∂∂ )()/(d 22 ( )+ dc u au uλ−∇ ∇ = ∂Ω ( )+n c u q u g∇ =r 1 1 1 1 2 2 1 2 1 1 2 2 2 2 + = + = h u h u r h u h u r 11 1 12 2 11 1 12 2 1 21 1 22 2 21 1 22 2 2 ( )+ ( )+ + = ( )+ ( )+ + = n c u n c u q u q u g n c u n c u q u q u g ∇ ∇ ∇ ∇ r r r r 11 1 12 2 1 11 1 12 2 11 1 12 2 1 11 21 1 22 2 21 1 22 2 2 12 + = ( )+ ( )+ + = + ( )+ ( )+ + = h u h u r n c u n c u q u q u g h n c u n c u q u q u g h µ µ ∇ ∇ ∇ ∇ r r r r ∂Ω nr ∂Ω 收稿日期：2005-09-20 作者简介：田兵（1984－），男，内蒙古包头人，学士，研究方向：数学物理方程。 万方数据 13 结构力学：平面应力 结构力学：平面应变 静电学 静磁学 交流电电磁学 直流电导电介质 热传导 扩散 确定问题类型后，可以在PDE Specification对话窗口中 填入c，a，f，d等系数（函数），这样就确定了待解的偏 微分方程。 1.2 确定边界条件 用函数assemb可以描述边界条件。 用pdetool提供的边界条件对话框，在对话框里填入g, h, q, r等边界条件。 1.3 确定偏微分方程所在Ω域的几何图形 平面上波的散射问题。 按照上面所说的解方程的过程，首先确定待解的偏微 分方程。 散射是介质对入射波的反射。假定介质是均匀的，那 么入射波在介质中传播的速度是一个常数c。波动方程通常 可以写成 其偏微分方程可以表示为 也可以写成 或者Helmholtz方程 k是波数，它与角频率ω，频率f，波长λ的关系是 可以用函数画出Ω域的几何图形。pdecirc：画圆； pdeellip：画椭圆；pderect：画矩形；pdepoly：画多边形。 也可以用鼠标在pdetool的画图窗口中直接画出Ω域的 几何图形。pdetool提供了类似于函数那样画圆、椭圆、矩 形、多边形的工具。 无论哪种画法，图形一经画出，pdetool就为这个图形 自动取名，并把代表图形的名字放入Set formula窗口，在 这个窗口，可以通过+、-图形的名字实现对图形的拓朴运 算，以便构造复杂的Ω域几何图形。 1.4 划分有限元 对Ω域进行有限元划分的函数有，initmesh：基本划 分；jigglemesh：采用jiggle方法划分；refinemesh：精细划 分。 在pdetool窗口中直接点击划分有限元的按钮划分有限 元，划分的方法与上面的函数相对应。 1.5 解方程 经过1.1~1.4步骤后就可以解方程。解方程的函数有， adaptmesh：解方程的通用函数；poicalc：矩形有限元快速 解椭圆型方程；poisolv：矩形有限元解椭圆型方程； parabolic：解抛物线型方程；hyperbolic：解双曲线型方 程。 在pdetool窗口中直接点击解方程的按钮即可解方程。 解方程所耗费的时间在于有限元划分的多少。 2 实例 我们现在设定边界条件。令入射光是在 =(cos(a)sin (a))方向上传播的平面波，其波动方程为 u是入射波v和反射波之和 u = v + r 边界条件比较简单 u = 0 因此 r = - v ( x , v ) 采用Sommerfeld散射条件。 当 趋向无穷，r可以表达为单向波方程 是反射方向，ξ是反射距离。 根据简谐波的条件，方程可以写成 ▽r=ikr 依据上面所说，可以在pdetool窗口中选择通用问题 （Generic Scalar），在PDE Specification对话框中填入c = 1，a = –k2 = –3600，f = 0。这样就确定了待解的偏微分方 程。 确定边界条件。 Dirichlet 条件：h=1， r=–exp(–i*60*x)。 确定偏微分方程所在域Ω的几何形状。 打开pdetool，画出被照射物体SQ1。SQ1是边长为0.1 的正方形，中心坐标是（0.5，0.5），将其旋转45度，令其 对角线分别平行于x，y坐标轴，x正方向上的顶点对着入射 波。我们考虑该顶点对入射波进行散射时的情况。我们的 计算域是圆形C1，其中心点与SQ1相同，半径0.5。这样， 偏微分方程所在域Ω的形状就是C1-SQ1，一个中心有着方 洞的圆形，类似于中国的古铜钱。 划分有限元。划分结果见图1。 图1 图2 解方程。解出的结果是复数，其实数部分如图2用颜色 表示，右边的标尺给出不同颜色所代表的数值。 pdetool可以多种方式画出计算结果，除了用颜色表示 外，还可以用等值线、箭头、立体图等表示。有限元结点 上的数值可以输出到文本。 -( , , ) ( , )e i tU x y t u x y ω= 2 2 2/ - =0U t c U∂ ∂ ∆ 2 2- - = 0u c uω ∆ 2- - = 0u k u∆ = / = 2 π / = 2 π /k c f cω λ a r ( - ) -( , , ) e ( , )ei k ax t iV x y t x yω ων= =r r ( , ) e ik a xx yν = r r x r / =0r t c rξ∂ ∂ + ∇r ξr ξr 万方数据 用MATLAB解偏微分方程 作者： 田兵， TIAN Bing 作者单位： 包头师范学院,学报编辑部,内蒙古,包头,014030 刊名： 阴山学刊（自然科学版） 英文刊名： YINSHAN ACADEMIC JOURNAL（NATURAL SCIENCE EDITION） 年，卷(期)： 2006,20(4) 被引用次数： 1次 引证文献(1条) 1.史策 热传导方程有限差分法的MATLAB实现[期刊论文]-咸阳师范学院学报 2009(4) 本文链接：http://d.g.wanfangdata.com.cn/Periodical_ysxk-z200604005.aspx