下载
加入VIP
  • 专属下载特权
  • 现金文档折扣购买
  • VIP免费专区
  • 千万文档免费下载

上传资料

关闭

关闭

关闭

封号提示

内容

首页 Matlab解微分方程

Matlab解微分方程.doc

Matlab解微分方程

阿泰斯特2013
2013-09-04 0人阅读 举报 0 0 暂无简介

简介:本文档为《Matlab解微分方程doc》,可适用于人文社科领域

第十六章偏微分方程的数值解法科学研究和工程技术中的许多问题可建立偏微分方程的数学模型。包含多个自变量的微分方程称为偏微分方程(partialdifferentialequation)简称PDE。偏微分方程问题其求解是十分困难的。除少数特殊情况外绝大多数情况均难以求出精确解。因此近似解法就显得更为重要。本章仅介绍求解各类典型偏微分方程定解问题的差分方法。几类偏微分方程的定解问题一个偏微分方程的表示通常如下:()式中是常数称为拟线性(quasilinear)数。通常存在种拟线性方程:双曲型(hyperbolic)方程:抛物线型(parabolic)方程:椭圆型(ellliptic)方程:。双曲型方程最简单形式为一阶双曲型方程:()物理中常见的一维振动与波动问题可用二阶波动方程:()描述它是双曲型方程的典型形式。方程的初值问题为:()边界条件一般有三类最简单的初边值问题为:()抛物型方程其最简单的形式为一维热传导方程:()方程可以有两种不同类型的定解问题:()初值问题:()()初边值问题:()其中为已知函数且满足连接条件:()边界条件为第一类边界条件。第二类和第三类边界条件为:()其中。当时为第二类边界条件否则称为第三类边界条件。椭圆型方程其最典型、最简单的形式是泊松(Poisson)方程()特别地当时即为拉普拉斯(Laplace)方程又称为调和方程:()Poisson方程的第一边值问题为:()其中为以为边界的有界区域为分段光滑曲线称为定解区域分别为上的已知连续函数。第二类和第三类边界条件可统一表示为:()其中为边界的外法线方向。当时为第二类边界条件时为第三类边界条件。差分方法的基本概念差分方法又称为有限差分方法或网格法是求偏微分方程定解问题的数值解中应用最广泛的方法之一。它的基本思想是:先对求解区域作网格剖分将自变量的连续变化区域用有限离散点(网格点)集代替将问题中出现的连续变量的函数用定义在网格点上离散变量的函数代替通过用网格点上函数的差商代替导数将含连续变量的偏微分方程定解问题化成只含有限个未知数的代数方程组(称为差分格式)。如果差分格式有解且当网格无限变小时其解收敛于原微分方程定解问题的解则差分格式的解就作为原问题的近似解(数值解)。因此用差分方法求偏微分方程定解问题一般需要解决以下问题:()选取网格()对微分方程及定解条件选择差分近似列出差分格式()求解差分格式()讨论差分格式解对于微分方程解的收敛性及误差估计。下面用一个简单的例子来说明用差分方法求解偏微分方程问题的一般过程及差分方法的基本概念。设有一阶双曲型方程初值问题。()选取网格:图差分示意图首先对定解区域作网格剖分最简单常用一种网格是用两族分别平行于轴与轴的等距直线:,()将分成许多小矩形区域。这些直线称为网格线其交点称为网格点也称为节点和分别称作方向和方向的步长。这种网格称为矩形网格。()对微分方程及定解条件选择差分近似列出差分格式。如果用向前差商表示一阶偏导数即:()其中。方程:()在节点处可表示为:()其中:。由于当足够小时在式中略去就得到一个与方程相近似的差分方程:()此处可看作是问题的解在节点处的近似值。同初值条件:()结合就得到求问题的数值解的差分格式。式:()称为差分方程的截断误差。如果一个差分方程的截断误差为则称差分方程对是阶精度对是阶精度的。显然截断误差的阶数越大差分方程对微分方程的逼近越好。若网格步长趋于时差分方程的截断误差也趋于则称差分方程与相应的微分方程是相容的。这是用差分方法求解偏微分方程问题的必要条件。如果当网格步长趋于时差分格式的解收敛到相应微分方程定解问题的解则称这种差分格式是收敛的。双曲型方程的差分解法一阶双曲型方程的差分格式考虑一阶双曲型方程的初值问题:()将平面剖分成矩形网格取方向步长为方向步长为网格线为:,为简便记:,以不同的差商近似偏导数可以得到方程的不同的差分近似()()()截断误差分别为与。结合离散化的初始条件可以得到几种简单的差分格式:()()()其中:。如果已知第层节点上的值按上面三种格式就可求出第层上的值。因此这三种格式都是显式格式。如果对采用向后差商采用向前差商则方程可化成:()相应的差分格式为:()此差分格式是一种隐式格式必须通过解方程组才能由第层节点上的值求出第层节点上的值。例对初值问题:其中:用差分格式求其数值解取。解:记由初始条件:按差分格式:计算公式为:。计算结果略。如果用差分格式:求解计算公式为:计算结果略。与准确解:比较知按前一个差分格式所求得的数值解不收敛到初值问题的解而后一个差分格式的解收敛到原问题的解。波动方程的差分格式对二阶波动方程:()如果令则方程可化成一阶线性双曲型方程组:()记则方程组可表成矩阵形式:()矩阵有两个不同的特征值故存在非奇异矩阵使得:()作变换方程组可化为:()方程组由二个独立的一阶双曲型方程联立而成。因此本节主要讨论一阶双曲型方程的差分解法。下面给出如下波动方程和边界条件的差分格式:()()将矩形划分成个小矩形长宽分别为:形成一个网格如图所示。图网格图可通过差分方程的方法计算近似值:在连续行内网格点的真实值为。求和的中心差分为:()()在每一行的网格间距是均匀的:而且。同时它在每一列也是均匀的。接下来将和去掉用()和()中的近似并按顺序代入()式可得到差分方程:()可用它来近似()式。为方便起见可将代入上式可得:()设行和的近似值已知可用上式求网格的行:()式中。根据上式左边的个已知值可得到近似值如图所示。图波动方程的空格样板用上式时必须注意如果计算的某个阶段带来的误差最终会越来越小则方法是稳定的。为了保证上式的稳定性必然使。还存在其他一些差分方程方法以称为隐格式法它们更难实现但对无限制。差分方法求解波动方程的MATLAB程序求解区间以()为边界条件的波动方程的差分方法程序。**********************************************************functionU=finedif(f,g,a,b,c,n,m)Inputf=u(x,)asastring'f'g=ut(x,)asastring'g'aandbrightendpointsof,aand,bctheconstantinthewaveequationnandmnumberofgridpointsover,aand,bOutputUsolutionmatrixanalogoustoTableInitializeparametersandUh=a(n)k=b(m)r=c*khr=r^r=r^s=r^s=*r^U=zeros(n,m)Computfirstandsecondrowsfori=:nU(i,)=feval(f,h*(i))U(i,)=s*feval(f,h*(i))k*feval(g,h*(i))r*(feval(f,h*i)feval(f,h*(i)))endComputeremainingrowsofUforj=:m,fori=:(n),U(i,j)=s*U(i,j)r*(U(i,j)U(i,j))U(i,j)endendU=U'*******************************************************************抛物型方程的差分解法以一维热传导方程:()为基本模型讨论适用于抛物型方程定解问题的几种差分格式。差分格式的建立首先对平面进行网格剖分。分别取为方向与方向的步长用两族平行直线将平面剖分成矩形网格节点为:。为简便记:。微分方程的差分近似在网格内点处对分别采用向前、向后及中心差商公式:()一维热传导方程可分别表示为:()()()由此得到一维热传导方程的不同差分近似:()()()上述差分方程所用到的节点各不相同。其截断误差分别为和。因此它们都与一维热传导方程相容。如果将式:中的用代替则可得到又一种差分近似:()差分方程用到四个节点。由Taylor公式容易得出:()故其的截断误差为。因而不是对任意的此差分方程都能逼近热传导方程:仅当时才成立。综上可知用不同的差商公式可以得到微分方程的不同的差分近似。构造差分格式的关键在于使其具有相容性、收敛性和稳定性。前面三个方程都具有相容性而此方程则要在一定条件下才具有相容性。初、边值条件的处理为用差分方法求解定解问题初值问题:()初边值问题:()还需对定解条件进行离散化。对初始条件及第一类边界条件可直接得到:式中:。对第二、三类边界条件:()需用差分近似。下面介绍两种较简单的处理方法。在左边界处用向前差商近似偏导数在右边界处用向后差商近似即:()则得边界条件的差分近似为:()其截断误差为。()用中心差商近似,即:()则得边界条件的差分近似为:()其截断误差为。误差的阶数提高了但出现定解区域外的节点和这就需要将解拓展到定解区域外。可以通过用内节点上的值插值求出和也可以假定热传导方程在边界上也成立将差分方程扩展到边界节点上由此消去和。几种常用的差分格式以热传导方程的初边值问题:()为例给出几种常用的差分格式。()古典显式格式令则:()可改写成:。将其与初始条件及第一类边界条件:结合我们得到求解此问题的一种差分格式:()由于第层上节点处的值已知由此即可算出在第一层上节点处的近似值。重复使用此式可以逐层计算出所有的因此此差分格式称为典显式格式。又因式中只出现相邻两个时间层的节点故此式是二层显式格式。()古典隐式格式将式:()整理并与初始条件及第一类边界条件式联立得差分格式如下:()其中:。虽然第层上的值仍为已知但不能由上式直接计算以上各层节点上的值必须通过解下列线性方程组:()式中。才能由计算故此差分格式称为古典隐式格式。此方程组是三对角方程组且系数矩阵严格对角占优故解存在唯一。()Richardson格式Richardson格式是将式:整理后与初始条件及第一类边界条件式联立。其计算公式为:)这种差分格式中所涉及的节点出现在三层上故为三层显式格式。Richardson格式是一种完全不稳定的差分格式因此它在实际计算中是不能采用的。()杜福特弗兰克尔(DoFortFrankel)格式DoFortFrankel格式也是三层显式格式它是由式:与初始条件及第一类边界条件式结合得到的。具体形式如下:()用这种格式求解时除了第层上的值由初值条件得到必须先用二层格式求出第层上的值然后再按上式逐层计算。()六点隐式格式对二阶中心差商公式:如果用在与处的二阶中心差商的平均值近似在处的值即:同时在点处的值也用中心差商:近似即:这样又得到热传导方程的一种差分近似:()其截断误差为将上式与初始条件及第一类边界条件式联立并整理得差分格式:()此格式涉及到六个节点它又是隐式格式故称为六点隐式格式。与古典隐式格式类似用六点格式由第层的值计算第层的值时需求解三对角方程组:()式中。此方程组的系数矩阵严格对角占优故仍可用追赶法求解。例用古典显式格式求初边值问题。的数值解取。解:这里:。由格式:可得到:将初值代入上式即可算出:将边界条件及上述结果代入又可求得:如此逐层计算得全部节点上的数值解为:前向差分法求解热传导方程的MATLAB程序设其中而且其中求解区间内的近似解。*****************************************************functionU=forwdif(f,c,c,a,b,c,n,m)Inputf=u(x,)asastring'f'c=u(,t)andc=u(a,t)aandbrightendpointsof,aand,bctheconstantintheheatequationnandmnumberofgridpointsover,aand,bOutputUsolutionmatrixanalogoustoTableInitializeparametersandUh=a(n)k=b(m)r=c^*kh^s=*rU=zeros(n,m)BoundaryconditionsU(,:m)=cU(n,:m)=cGeneratefirstrowU(:n,)=feval(f,h:h:(n)*h)'GenerateremainingrowsofUforj=:mfori=:nU(i,j)=s*U(i,j)r*(U(i,j)U(i,j))endendU=U'*****************************************************CrankNicholson求解热传导方程的MATLAB程序设其中而且其中求解区间内的近似解。*****************************************************functionU=crnich(f,c,c,a,b,c,n,m)Inputf=u(x,)asastring'f'c=u(,t)andc=u(a,t)aandbrightendpointsof,aand,bctheconstantintheheatequationnandmnumberofgridpointsover,aand,bOutputUsolutionmatrixanalogoustoTableInitializeparametersandUh=a(n)k=b(m)r=c^*kh^s=rs=rU=zeros(n,m)BoundaryconditionsU(,:m)=cU(n,:m)=cGeneratefirstrowU(:n,)=feval(f,h:h:(n)*h)'FormthediagonalandoffdiagonalelemntsofAandtheconstantvectorBandsolvetridiagonalsystemAX=BVd(,:n)=s*ones(,n)Vd()=Vd(n)=Va=ones(,n)Va(n)=Vc=ones(,n)Vc()=Vb()=cVb(n)=cforj=:mfori=:nVb(i)=U(i,j)U(i,j)s*U(i,j)endX=trisys(Va,Vd,Vc,Vb)U(:n,j)=X'endU=U'*****************************************************椭圆型方程第一边值问题的差分解法本节以Poisson方程为基本模型讨论第一边值问题的差分方法。差分格式的建立考虑Poisson方程的第一边值问题:()图取分别为方向和方向的步长如图所示以两族平行线:将定解区域剖分成矩形网格。节点的全体记为:。定解区域内部的节点称为内点记内点集为。边界与网格线的交点称为边界点边界点全体记为。与节点沿方向或方向只差一个步长的点和称为节点的相邻节点。如果一个内点的四个相邻节点均属于称为正则内点正内点的全体记为至少有一个相邻节点不属于的内点称为非正则内点非正则内点的全体记为。问题是要求出第一边值问题在全体内点上的数值解。为简便记:。对正则内点由二阶中心差商公式:()Poisson方程在点处可表示为:()其中:()为其截断误差表示式略去即得与方程相近似的差分方程:()式中方程的个数等于正则内点的个数而未知数则除了包含正则内点处解的近似值外还包含一些非正则内点处的近似值因而方程个数少于未知数个数。在非正则内点处Poisson方程的差分近似不能按上式给出需要利用边界条件得到。边界条件的处理可以有各种方案下面介绍较简单的两种。()直接转移。用最接近非正则内点的边界点上的值作为该点上值的近似这就是边界条件的直接转移。例如点为非正则内点其最接近的边界点为点则有:()上式可以看作是用零次插值得到非正则内点处的近似值容易求出其截断误差为。将上式代入方程个数即与未知数个数相等。()线性插值。这种方案是通过用同一条网格线上与点相邻的边界点与内点作线性插值得到非正则内点处值的近似。由点与的线性插值确定的近似值得:()其中其截断误差为。将其与方程相近似的差分方程联立得到方程个数与未知数个数相等的方程组求解此方程组可得Poisson方程第一边值问题的数值解。上面所给出的差分格式称为五点菱形格式()实际计算时经常取此时五点菱形格式可化为:()简记为:()其中:。例用五点菱形格式求解拉普拉斯(Laplace)方程第一边值问题。其中:。取。解:网格中有四个内点均为正则内点。由五点菱形格式得方程组:代入边界条件:其解为:当时对利用点构造的差分格式称为五点矩形格式简记为:()其中:其截断误差为:()五点菱形格式与矩形格式的截断误差均为称它们具有二阶精度。如果用更多的点构造差分格式其截断误差的阶数可以提高如利用菱形格式及矩形格式所涉及的所有节点构造出的九点格式就是具有四阶精度的差分格式。用Dirichlet法求解Laplace方程的MATLAB程序求解区间内的近似解。而且满足条件:其中而且其中而且存在整数和使得:。*****************************************************functionU=dirich(f,f,f,f,a,b,h,tol,max)Inputf,f,f,fareboundaryfunctionsinputasstringsaandbrightendpointsof,aand,bhstepsizetolisthetoleranceOutputUsolutionmatrixanalogoustoTableInitializeparametersandUn=fix(ah)m=fix(bh)ave=(a*(feval(f,)feval(f,))b*(feval(f,)feval(f,)))(*a*b)U=ave*ones(n,m)BoundaryconditionsU(,:m)=feval(f,:h:(m)*h)'U(n,:m)=feval(f,:h:(m)*h)'U(:n,)=feval(f,:h:(n)*h)U(:n,m)=feval(f,:h:(n)*h)U(,)=(U(,)U(,))U(,m)=(U(,m)U(,m))U(n,)=(U(n,)U(n,))U(n,m)=(U(n,m)U(n,m))SORparameterw=(sqrt((cos(pi(n))cos(pi(m)))^))Refineapproximationsandsweepoperatorthroughoutthegriderr=cnt=while((err>tol)(cnt<=max))err=forj=:mfori=:nrelx=w*(U(i,j)U(i,j)U(i,j)U(i,j)*U(i,j))U(i,j)=U(i,j)relxif(err<=abs(relx))err=abs(relx)endendendcnt=cntendU=flipud(U')*****************************************************unknownunknownunknownunknownunknownunknownunknownunknownunknownunknownunknownunknownunknownunknownunknownunknownunknownunknownunknownunknownunknownunknownunknownunknownunknownunknownunknownunknownunknownunknownunknownunknownunknownunknownunknownunknownunknownunknownunknownunknownunknownunknownunknownunknownunknownunknownunknownunknownunknownunknownunknownunknownunknownunknownunknownunknownunknownunknownunknownunknownunknownunknownunknownunknownunknownunknownunknownunknownunknownunknownunknownunknownunknownunknownunknownunknownunknownunknownunknownunknownunknownunknownunknownunknownunknownunknownunknownunknownunknownunknownunknownunknownunknownunknownunknownunknownunknownunknownunknownunknownunknownunknownunknownunknownunknownunknownunknownunknownunknownunknownunknownunknownunknownunknownunknownunknownunknownunknownunknownunknownunknownunknownunknownunknownunknownunknownunknownunknownunknownunknownunknownunknownunknownunknownunknownunknownunknownunknownunknownunknownunknownunknownunknownunknownunknownunknownunknownunknownunknownunknownunknownunknownunknownunknownunknownunknownunknownunknownunknownunknownunknownunknownunknownunknownunknownunknownunknownunknownunknownunknownunknownunknownunknownunknownunknownunknownunknownunknownunknownunknownunknownunknownunknownunknownunknownunknownunknownunknownunknownunknownunknownunknownunknownunknownunknownunknownunknownunknownunknownunknownunknownunknownunknownunknownunknownunknownunknownunknownunknownunknownunknownunknownunknownunknownunknownunknownunknownunknownunknownunknownunknownunknownunknownunknownunknownunknownunknownunknownunknownunknownunknownunknownunknownunknownunknownunknownunknownunknownunknownunknownunknownunknownunknownunknownunknownunknownunknownunknownunknownunknownunknownunknownunknownunknownunknownunknownunknownunknownunknownunknownunknownunknownunknownunknownunknownunknownunknownunknownunknownunknownunknownunknownunknownunknownunknownunknownunknownunknownunknown

用户评价(0)

关闭

新课改视野下建构高中语文教学实验成果报告(32KB)

抱歉,积分不足下载失败,请稍后再试!

提示

试读已结束,如需要继续阅读或者下载,敬请购买!

文档小程序码

使用微信“扫一扫”扫码寻找文档

1

打开微信

2

扫描小程序码

3

发布寻找信息

4

等待寻找结果

我知道了
评分:

/27

Matlab解微分方程

VIP

在线
客服

免费
邮箱

爱问共享资料服务号

扫描关注领取更多福利