偏微分方程—matlab

基础知识偏微分方程的定解问题各种物理性质的定常（即不随时间变化）过程，都可用椭圆型方程来描述。其最典型、最简单的形式是泊松(Poisson)方程),(2222yxfyuxuu（1）特别地，当f(x,y)≡0时，即为拉普拉斯(Laplace)方程，又称为调和方程02222yuxuu（2）带有稳定热源或内部无热源的稳定温度场的温度分布，不可压缩流体的稳定无旋流动及静电场的电势等均满足这类方程。Poisson方程的第一边值问题为),(),(),(),(),(2222yxyxuyxyxfyuxuuyx（3）其中Ω为以Γ为边界的有界区域，Γ为分段光滑曲线，ΩUΓ称为定解区域，f(x,y),ϕ(x,y)分别为Ω,Γ上的已知连续函数。第二类和第三类边界条件可统一表示成)0(0),(aunuyx（4）其中n为边界Γ的外法线方向。当α=0时为第二类边界条件，α≠0时为第三类边界条件。在研究热传导过程，气体扩散现象及电磁场的传播等随时间变化的非定常物理问题时，常常会遇到抛物型方程。其最简单的形式为一维热传导方程)0(022axuatu（5）方程（5）可以有两种不同类型的定解问题：初值问题（也称为Cauchy问题）xxxuxtxuatu)()0,(,0022（6）初边值问题lxtgtlutgtuxxulxTtxuatu0),(),(),(),0()()0,(0,002122（7）其中ϕ)(),(),(21xgxgx为已知函数，且满足连接条件)0()(),0()0(21glg问题（7）中的边界条件)(),(),(),0(21tgtlutgtu称为第一类界条件。第二类和第三类边界条件为TttgutxuTttgutxulxx0),()(0),()(22101（8）其中0,021。当021时，为第二类边界条件，否则称为第三类边界条件。双曲型方程的最简单形式为一阶双曲型方程0xuatu（9）物理中常见的一维振动与波动问题可用二阶波动方程2222xuatu（10）描述，它是双曲型方程的典型形式。方程（10）的初值问题为xxtuxxxuxtxuatut)()()0,(,002222（11）边界条件一般也有三类，最简单的初边值问题为Tttgtlutgtulxxtuxxulxtxuatut0)(),(),(),0(0)()()0,(0,02102222如果偏微分方程定解问题的解存在，唯一且连续依赖于定解数据（即出现在方程和定解条件中的已知函数），则此定解问题是适定的。可以证明，上面所举各种定解问题都是适定的。§2偏微分方程的差分解法差分方法又称为有限差分方法或网格法，是求偏微分方程定解问题的数值解中应用最广泛的方法之一。它的基本思想是：先对求解区域作网格剖分，将自变量的连续变化区域用有限离散点（网格点）集代替；将问题中出现的连续变量的函数用定义在网格点上离散变量的函数代替；通过用网格点上函数的差商代替导数，将含连续变量的偏微分方程定解问题化成只含有限个未知数的代数方程组（称为差分格式）。如果差分格式有解，且当网格无限变小时其解收敛于原微分方程定解问题的解，则差分格式的解就作为原问题的近似解（数值解）。因此，用差分方法求偏微分方程定解问题一般需要解决以下问题：（i）选取网格；（ii）对微分方程及定解条件选择差分近似，列出差分格式；（iii）求解差分格式；（iv）讨论差分格式解对于微分方程解的收敛性及误差估计。下面我们只对偏微分方程的差分解法作一简要的介绍。2.1椭圆型方程第一边值问题的差分解法以Poisson方程（1）为基本模型讨论第一边值问题的差分方法。考虑Poisson方程的第一边值问题（3）),(),(),(),(),(2222yxyxuyxyxfyuxuyx取h,τ分别为x方向和y方向的步长，以两族平行线jyykhxxjk,),2,1,0,(jk将定解区域剖分成矩形网格。节点的全体记为,,|),{(jykhxyxRjkkkji,为整数}。定解区域内部的节点称为内点，记内点集R为h。边界Γ与网格线的交点称为边界点，边界点全体记为Γhτ。与节点),(jkyx沿x方向或y方向只差一个步长的点),(1jkyx和),(1jkyx称为节点),(jkyx的相邻节点。如果一个内点的四个相邻节点均属于ΩUΓ，称为正则内点，正则内点的全体记为Ω(1)，至少有一个相邻节点不属于ΩUΓ的内点称为非正则内点，非正则内点的全体记为Ω(2)。我们的问题是要求出问题（3）在全体内点上的数值解。为简便记，记),(),,(),(),(),(,jkjkjkjkyxffyxujkuyxjk。对正则内点)1(),(jk，由二阶中心差商公式)(),1(),(2),1(22),(22hOhjkujkujkuxujk)()1,(),(2)1,(22),(22OjkujkujkuyujkPoisson方程（1）在点),(jk处可表示为)(2222,21,,1,2,1,,1hOfuuuhuuujkjkjkjkjkjkjk（12）在式（12）中略去)(22hO，即得与方程（1）相近似的差分方程jkjkjkjkjkjkjkfuuuhuuu,21,,1,2,1,,122（13）式（13）中方程的个数等于正则内点的个数，而未知数jku,,则除了包含正则内点处解u的近似值，还包含一些非正则内点处u的近似值，因而方程个数少于未知数个数。在非正则内点处Poisson方程的差分近似不能按式（13）给出，需要利用边界条件得到。边界条件的处理可以有各种方案，下面介绍较简单的两种。(i)直接转移(ii)线性插值由式（13）所给出的差分格式称为五点菱形格式，实际计算时经常取h=τ，此时五点菱形格式可化为jkjkjkjkjkjkfuuuuuh,,1,1,,1,1241（14）简记为21hjkjkfu,,（15）其中jkjkjkjkjkjkuuuuuu,1,1,,1,1,4。求解差分方程组最常用的方法是同步迭代法，同步迭代法是最简单的迭代方式。除边界节点外，区域内节点的初始值是任意取定的。例1用五点菱形格式求解Laplace方程第一边值问题22),(2222)1(lg),(),(0yxyxuyxyuxuyx其中}1,0|),{(yxyx。取31h。当h时，利用点(k,j),(k±1,j.1),(k±1,j+1)构造的差分格式jkjkjkjkjkjkfuuuuuh,,1,11,11,11,12421（16）称为五点矩形格式，简记为221hjkjkfu,,（17）其中jkjkjkjkjkjkuuuuuu,1,11,11,11,1,4。2.2抛物型方程的差分解法以一维热传导方程（5）)0(022axuatu为基本模型讨论适用于抛物型方程定解问题的几种差分格式。首先对xt平面进行网格剖分。分别取h,τ为x方向与t方向的步长，用两族平行直线khxxk(k=0,±1,±2,…)，jttjk(j=0,1,2,…)，将xt平面剖分成矩形网格，节点为),(jkyx(k=0,±1,±2,…,j=0,1,2,…)。为简便起见，记),,(),(jkyxjk),,(),(jkyxujku),(kkx),(11jjtgg),(22jjtgg),(11jjt)(22jjt。2.2.1微分方程的差分近似在网格内点(k,j)处，对tu分别采用向前、向后及中心差商公式，对22xu采用二阶中心差商公式，一维热传导方程（5）可分别表示为)(),1(),(2),1(2)1,()1,()(),1(),(2),1()1,(),()(),1(),(2),1(),()1,(222222hOhjkujkujkuajkujkuhOhjkujkujkuajkujkuhOhjkujkujkuajkujku由此得到一维热传导方程的不同的差分近似022,1,,1,1,huuuauujkjkjkjkjk（18）022,1,,11,,huuuauujkjkjkjkjk（19）0222,1,,11,1,huuuauujkjkjkjkjk（20）2.2.2初、边值条件的处理为用差分方程求解定解问题（6），（7）等，还需对定解条件进行离散化。对初始条件及第一类边界条件，可直接得到)1,0,1,0(0,0,nkkxuukkk或（21）jjjnjjjgtluugtuu2,1,0),(),0()1,1,0(mj（22）其中Tmhln,。对第二、三类边界条件则需用差商近似。下面介绍两种较简单的处理方法。（i）在左边界(x=0)处用向前差商近似偏导数xu，在右边界(lx)处用向后差商近似偏导数xu，即即得边界条件（8）的差分近似为（ii）用中心差商近似xu，即则得边界条件的差分近似为这样处理边界条件，误差的阶数提高了，但式（24）中出现定解区域外的节点(-1,j)和(n+1,j)，这就需要将解拓展到定解区域外。可以通过用内节点上的u值插值求出和，也可以假定热传导方程（5）在边界上也成立，将差分方程扩展到边界节点上，由此消去和。2.2.3几种常用的差分格式下面我们以热传导方程的初边值问题（7）为例给出几种常用的差分格式。(i)古典显式格式为便于计算，令，式（18）改写成以下形式将式（18）与（21），（22）结合，我们得到求解问题（7）的一种差分格式由于第0层(j=0)上节点处的u值已知，由式（25）即可算出u在第一层(j=1)上节点处的近似值。重复使用式（25），可以逐层计算出各层节点的近似值。（ii）古典隐式格式将（19）整理并与式（21），（22）联立，得差分格式如下其中。虽然第0层上的u值仍为已知，但不能由式（30）直接计算以上各层节点上的值故差分格式（26）称为古典隐式格式。（iii）杜福特—弗兰克尔（DoFort—Frankel）格式DoFort—Frankel格式是三层显式格式，它是由式（24）与（25），（26）结合得到的。具体形式如下：用这种格式求解时，除了第0层上的值由初值条件（21）得到，必须先用二层格式求出第1层上的值，然后再按格式（27）逐层计算。2.3双曲型方程的差分解法对二阶波动方程（10）如果令，则方程（10）可化成一阶线性双曲型方程组记，则方程组（28）可表成矩阵形式矩阵A有两个不同的特征值λ=±a，故存在非奇异矩阵P，使得作变换，方程组（29）可化成方程组（30）由两个独立的一阶双曲型方程联立而成。因此下面主要讨论一阶双曲型方程的差分解法。考虑一阶双曲型方程的初值问题与抛物型方程的讨论类似，仍将xt平面剖分成矩形网格。取x方向步长为h，t方向步长为τ，网格线为。为简便起见，记。以不同的差商近似偏导数，可以得到方程（9）的不同的差分近似结合离散化的初始条件，可以得到几种简单的差分格式。§3一维状态空间的偏微分方程的MATLAB解法3.1工具箱命令介绍MATLAB提供了一个指令pdepe，用以解以下的PDE方程式其中时间介于之间，而