《偏微分方程概述及运用matlab求解偏微分方程常见问题》Word格式.docx
《《偏微分方程概述及运用matlab求解偏微分方程常见问题》Word格式.docx》由会员分享,可在线阅读,更多相关《《偏微分方程概述及运用matlab求解偏微分方程常见问题》Word格式.docx(18页珍藏版)》请在冰豆网上搜索。
在我国,偏微分方程的研究起步较晚。
但解放后,在党和国家的大力号召和积极支持下,我国偏微分方程的研究工作发展比较迅速,涌现出一批在这一领域中做出杰出工作的数学家,如谷超豪院士、李大潜院士等,并在一些研究方向上达到了国际先进水平。
但总体来说,偏微分方程的研究队伍的组织和水平、研究工作的广度和深度与世界先进水平相比还有很大的差距。
因此,我们必须继续努力,大力加强应用偏微分方程的研究,逐步缩小与世界先进水平的差距
二,偏微分方程的内容
偏微分方程是什么样的?
它包括哪些内容?
这里我们可从一个例子的研究加以介绍。
弦振动是一种机械运动,当然机械运动的基本定律是质点力学的F=ma,但是弦并不是质点,所以质点力学的定律并不适用在弦振动的研究上。
然而,如果我们把弦细细地分成若干个极小极小的小段,每一小段抽象地看作是一个质点,这样我们就可以应用质点力学的基本定律了。
弦是指又细又长的弹性物质,比如弦乐器所用的弦就是细长的、柔软的、带有弹性的。
演奏的时候,弦总是绷紧着具有一种张力,这种张力大于弦的重量几万倍。
当演奏的人用薄片拨动或者用弓在弦上拉动,虽然只因其所接触的一段弦振动,但是由于张力的作用,传播到使整个弦振动起来。
用微分的方法分析可得到弦上一点的位移是这一点所在的位置和时间为自变量的偏微分方程。
偏方程又很多种类型,一般包括椭圆型偏微分方程、抛物型偏微分方程、双曲型偏微分方程。
上述的例子是弦振动方程,它属于数学物理方程中的波动方程,也就是双曲型偏微分方程。
偏微分方程的解一般有无穷多个,但是解决具体的物理问题的时候,必须从中选取所需要的解,因此,还必须知道附加条件。
因为偏微分方程是同一类现象的共同规律的表示式,仅仅知道这种共同规律还不足以掌握和了解具体问题的特殊性,所以就物理现象来说,各个具体问题的特殊性就在于研究对象所处的特定条件,就是初始条件和边界条件。
拿上面所举的弦振动的例子来说,对于同样的弦的弦乐器,如果一种是以薄片拨动弦,另一种是以弓在弦上拉动,那么它们发出的声音是不同的。
原因就是由于“拨动”或“拉动”的那个“初始”时刻的振动情况不同,因此产生后来的振动情况也就不同。
天文学中也有类似情况,如果要通过计算预言天体的运动,必须要知道这些天体的质量,同时除了牛顿定律的一般公式外,还必须知道我们所研究的天体系统的初始状态,就是在某个起始时间,这些天体的分布以及它们的速度。
在解决任何数学物理方程的时候,总会有类似的附加条件。
就弦振动来说,弦振动方程只表示弦的内点的力学规律,对弦的端点就不成立,所以在弦的两端必须给出边界条件,也就是考虑研究对象所处的边界上的物理状况。
边界条件也叫做边值问题。
当然,客观实际中也还是有“没有初始条件的问题”,如定场问题(静电场、稳定浓度分布、稳定温度分布等),也有“没有边界条件的问题”,如着重研究不靠近两端的那段弦,就抽象的成为无边界的弦了。
在数学上,初始条件和边界条件叫做定解条件。
偏微分方程本身是表达同一类物理现象的共性,是作为解决问题的依据;
定解条件却反映出具体问题的个性,它提出了问题的具体情况。
方程和定解条件合而为一体,就叫做定解问题。
求偏微分方程的定解问题可以先求出它的通解,然后再用定解条件确定出函数。
但是一般来说,在实际中通解是不容易求出的,用定解条件确定函数更是比较困难的。
偏微分方程的解法还可以用分离系数法,也叫做傅立叶级数;
还可以用分离变数法,也叫做傅立叶变换或傅立叶积分。
分离系数法可以求解有界空间中的定解问题,分离变数法可以求解无界空间的定解问题;
也可以用拉普拉斯变换法去求解一维空间的数学物理方程的定解。
对方程实行拉普拉斯变换可以转化成常微分方程,而且初始条件也一并考虑到,解出常微分方程后进行反演就可以了。
应该指出,偏微分方程的定解虽然有以上各种解法,但是我们不能忽视由于某些原因有许多定解问题是不能严格解出的,只可以用近似方法求出满足实际需要的近似程度的近似解。
常用的方法有变分法和有限差分法。
变分法是把定解问题转化成变分问题,再求变分问题的近似解;
有限差分法是把定解问题转化成代数方程,然后用计算机进行计算;
还有一种更有意义的模拟法,它用另一个物理的问题实验研究来代替所研究某个物理问题的定解。
虽然物理现象本质不同,但是抽象地表示在数学上是同一个定解问题,如研究某个不规则形状的物体里的稳定温度分布问题,在数学上是拉普拉斯方程的边值问题,由于求解比较困难,可作相应的静电场或稳恒电流场实验研究,测定场中各处的电势,从而也解决了所研究的稳定温度场中的温度分布问题。
随着物理科学所研究的现象在广度和深度两方面的扩展,偏微分方程的应用范围更广泛。
从数学自身的角度看,偏微分方程的求解促使数学在函数论、变分法、级数展开、常微分方程、代数、微分几何等各方面进行发展。
从这个角度说,偏微分方程变成了数学的中心。
三,用matlab解偏微分方程
解偏微分方程不是一件轻松的事,但是偏微分方程在自然科学和工程领域中应用很广,因此,我们可以运用matlab这个软件来解决一些常见的偏微分方程。
(一)Matlab偏微分方程工具箱简介。
1,概述。
本文只给出该工具箱的函数列表
2,偏微分方程算法函数列表。
adaptmesh生成自适应网络及偏微分方程的解
assemb生成边界质量和刚度矩阵
assema生成积分区域上质量和刚度矩阵
assempde组成偏微分方程的刚度矩阵及右边
hyperbolic求解双曲线型偏微分方程
parabolic求解抛物线型偏微分方程
pdeeig求解特征型偏微分方程
pdenonlin求解非线性型微分方程
poisolv利用矩阵格式快速求解泊松方程
3,图形界面函数。
pdecirc画圆
pdeellip画椭圆
pdemdlcv转化为版本1.0式的*.m文件
pdepoly画多边形
pderect画矩形
pdetool偏微分方程工具箱的图形用户界面
4,几何处理函数。
csgchk检查几何矩阵的有效性
csgdel删除接近边界的小区
decsg将固定的几何区域分解为最小区域
initmesh产生最初的三角形网络
jigglemesh微调区域内的三角形网络
poimesh在矩形区域上产生规则的网络
refinemesh细化三角形网络
wbound写一个边界描述文件
wgeom写一个几何描述文件
pdecont画轮廓图
pdemesh画偏微分方程的三角形网络
pdeplot画偏微分方程的三角形网络
pdesurf画表面图命令
5,通用函数。
pdetriq三角形单元的品性度量
poiasma边界点对快速求解泊松方程的“贡献”矩阵
poicalc规范化的矩阵格式的点索引
(二)Matlab偏微分方程工具箱应用。
可以用词工具箱求解如椭圆方程,双曲线方程,特征值方程,抛物线方程。
椭圆型偏微分方程
椭圆型偏微分方程的一般形式为
其中:
若
,
为
的梯度,则其定义为
散度
的定义为
这样,
可以更明确地表示为
为常数,则进一步化简为
其中,
又称为Laplace算子。
这样椭圆型偏微分方程可以简单地写为
抛物型偏微分方程
抛物型偏微分方程的一般形式为
根据上面叙述,若
为常数,则该方程可以更简单地写为
双曲型偏微分方程
双曲型偏微分方程的一般形式为
为常数,则可以将该方程简化为
三类方程的直接的区别在于
对
的导数的阶次。
若对
没有求导,可以理解为其值为常数,故称为椭圆型的。
若取
对时间
的一阶导数,则与
的二阶导数直接构成了抛物线关系,故称为抛物型偏微分方程。
的二阶导数,称其为双曲型偏微分方程。
特征值型偏微分方程
特征值型偏微分方程为
对常数
该方程可以化简为
该方程是椭圆型偏微分方程的特例。
pdetool的使用:
在matlab命令窗口中键入pdetool窗口打开进入工作状态,pdetool提供两种解方程的方法,一种是通过函数,利用函数可以编程也可以用命令行的方式解方程,两一种是对pdetool窗口进行交互操作。
一般来说,用函数解方程比较繁琐,但是比较灵活:
通过窗口交互操作比较简单。
解方程的全部过程以及结果都可以输出保存为文本文件。
限于文本的篇幅,我们主要介绍交互操作偏微分方程的方法。
1.确定待解的偏微分方程。
利用函数assempde可以对待解的偏微分方程加以描述。
在交互操作中,为了方便用户,pdetool把常见问题归结为及各类型,可以再pdetool窗口的工具栏上找到选择类型的弹出菜单,这些类型如下:
通用问题
通用系统(二维的偏微分方程组)
结构力学:
平面应力
平面应变
静电学
静磁学
交流电电磁学
直流电导电介质
热传导
扩散
确定问题类型后,可以再PDESpecification对话框填入c,a,f,d等系数,这样就确定了待解的偏微分方程。
2.确定边界条件
用函数assemb可以描述边界条件。
用pdetool提供的边界条件对话框,在对话框里填入g,h,q,r等边界条件。
3.确定偏微分方程所在
域的几何图形
平面上波得散射问题。
按照上面所说的解方程的过程,首先确定带解的偏微分方程。
散射是介质对入射波的反射。
假定介质是均匀的,那么入射波在介质中传播的速度是一个常数c。
可以用函数画出
域的几何图形。
Pdecirc:
画圆:
pdeellip:
画椭圆:
pderect:
画矩形:
pdepoly:
画多边形。
无论哪种画法,图形一经画出,pdetool就为这个图形自动取名,并把代表图形的名字放入Setformula窗口,在这个窗口,可以通过+,-图形的名字现在对图形的拓扑运算,以便构造复杂的
域几何图形
4.划分有限元
域进行有限元的划分函数有,initmesh:
基本划分;
jigglemesh:
采用jiggle方法划分;
refinemesh:
精细划分。
在pdetool窗口中直接点击划分有限元的按钮划分有限元,划分的方法与上面的函数想对应。
5.解方程
经过1.1—1.4步骤后就可以解方程。
解方程的函数有,adaptunesh:
解方程的通用函数;
poicalc:
矩形有限元快速解椭圆形方程;
poisolv:
矩形有限元解椭圆形方程;
parabolic:
解抛物线型方程:
hyperbolic:
解双曲线型方程。
在pdetool窗口中直接点击解方程的按钮即可解方程,解方程所耗费的时间在于有限元划分的多少。
(三)实例
求解偏微分方程的边值问题
一、MATLAB支持的偏微分方程类型
考虑平面有界区域D上的二阶椭圆型PDE边值问题:
其中
未知函数为
。
它的边界条件分为三类:
(1)Direchlet条件:
(2)Neumann条件:
(3)混合边界条件:
在边界
上部分为Direchlet条件,另外部分为Neumann条件。
是定义在边界
的已知函数,另外
也可以是一个2*2的函数矩阵,
是沿边界的外法线的单位向量。
在使用pdetool时要向它提供这些已知参数。
二、例题
例题1用pdetool求解
解:
首先在MATLAB的工作命令行中键入pdetool,按回牟键确定,于是出现PDEToolbox窗口,选GenenicScalar模式.
(l)画区域圆
单击椭圆工具按钮,大致在(0,0)位置单击鼠标右键,拖拉鼠标到适当位置松开。
为了保证所绘制的圆是标准的单位园,在所绘园上双击,打开ObjectDialog对话框,精确地输入圆心坐标X-center为0、Y-center为0及半径Radius为l,然后单击OK按钮,这样单位画已画好.
(2)设置边界条件
单击工具边界模式按钮,图形边界变红,逐段双击边界,打开Boundarycondition对话框.输入边界条件.对于同一类型的边界,可以按Shift键,将多个边界同时选择,统一设边界条件.本题选择Dirichlet条件,输入h为1,r为0。
,然后单击OK按钮.也可以单击Boundary菜单中SpocifyBoundaryCondition…选项,打开BoundaryCondition对话框输入边界条件.
(3)设置方程
单击偏微分方程按钮,打开PDESpecification对话框,选择方程类型·
本题选Ellintic(椭圆型),输入c为1,a为O,f为1,然后单击OK按钮.
(4)网格剖分
单击网格工具,或者单击Mesh菜单中InitializeMesh项,可进行初始网格剖分.这时在PDEToolbox窗口下方的状态栏内显示出初始网格的节点数和三角形单元数.本题节点数为144个,三角形单元数为254个(图?
?
)。
如果要细化网格,单击细化工具,或者单击Mesh菜单中RefineMesh选项,节点数成为541个,三角形单元数为1016个。
(5)解方程
单击解方程工具,或者单击Solve菜单中SolvePDE选项,可求得方程数值解并用彩色图形显示。
单击作图工具,或者单击Plot菜单中Parameter…选项,出现Plotselection对话框.从中选择于Height(3-Dplot),然后单击Plot按钮,方程的图形解如图?
所示。
除了作定解问题解u的图形外,也可以作
等图形·
(6)输出网格节点的编号、单元编号以及节点坐标
单击Mesh菜单中ShowNodeLabels选项,再单击网格工具,即可显示节点编号(图?
?
)。
若要输出节点坐标,只需单击Mesh菜单中ExportMesh…选项,这时打开的Export对话框中的默认值为pet,这里p、e、t分别表示point(点)、edges(边)、triangles(三角形)数据变量,单击OK按钮,然后在MATLAB命令行键入p,即可以显示按节点编号排列的坐标;
键入e再回车则显示边界数据矩阵(7维数组);
键入t按回车则显示三角形单元数据矩阵(4维数组)。
点、边、单元的部分输出为:
p=
Columns1through11
-1.00000.00001.00000.0000-0.70710.70710.7071-0.7071-0.9808-0.9239-0.8315
-0.0000-1.000001.0000-0.7071-0.70710.70710.7071-0.1951-0.3827-0.5556
e=
1.00009.000010.000011.00005.000012.000013.000014.00002.000015.000016.0000
9.000010.000011.00005.000012.000013.000014.00002.000015.000016.000017.0000
00.12500.25000.37500.50000.62500.75000.875000.12500.2500
0.12500.25000.37500.50000.62500.75000.87501.00000.12500.25000.3750
1.00001.00001.00001.00001.00001.00001.00001.00002.00002.00002.0000
1.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.0000
00000000000
t=
Columns1through18
3214202629172310011668919945121397
123486728532191051213142
899781988492991279489119119951181189070126
111111111111111111
(7)输出近似解数值
单击Solve菜单中ExportSolution。
选项,在Export对话框中输入u再单击OK按钮,再在MATLAB命令行中输入u并回车,就会显示按节点编号排列的解u的数值。
(8)近似解和准确解的比较
方程的准确解为:
为了与准确解比较,单击Plot菜单中Parameters…选项,打开PlotSelection对话框,在Height(3-Dplot)行的Property下拉框中选UserEntry,并且输入
u-(1-x.^2-y.^2)/4,单击Plot按钮,就可以看到误差曲面,其数量级为
参考文献
1,《偏微分方程的MATLAB解法》(作者:
陆君安)
2,《用MATLAB解偏微分方程》(阴山学刊:
作者:
田兵)
3,XX百科