大连理工大学-高等数值分析-抛物型方程有限差分法

抛物型方程有限差分法抛物方程差分法的构造在空间方向上与椭圆方程类似，在时间方向上用一阶差商代替代替一阶微商。然后在时间方向上逐层求解。特别当空间维数较高时，可以使用局部一维格式大大降低计算量。1.简单差分法考虑一维模型热传导方程(1.1))(22xfxuatu，Tt0其中a为常数。)(xf是给定的连续函数。（1.1）的定解问题分两类：第一，初值问题（Cauchy问题）：求足够光滑的函数txu,，满足方程（1.1）和初始条件：（1.2）xxu0,,x第二，初边值问题（也称混合问题）：求足够光滑的函数txu,，满足方程（1.1）和初始条件：13.1xxu0,,lxl及边值条件23.10,,0tlutu，Tt0假定xf和x在相应的区域光滑，并且于0,0，0,l两点满足相容条件,则上述问题有唯一的充分光滑的解。现在考虑边值问题（1.1），（1.3）的差分逼近取Nlh为空间步长，MT为时间步长，其中N，M是自然数，jhxxj,Nj,,1,0；kyyk,Mk,,1,0将矩形域GTtlx0;0分割成矩形网格。其中jiyx,表示网格节点；hG表示网格内点（位于开矩形G中的网格节点）的集合；hG表示位于闭矩形G中的网格节点的集合；h表示hG-hG网格边界点的集合。kju表示定义在网点kitx,处的待求近似解，Nj0，Mk0。注意到在节点kitx,处的微商和差商之间的下列关系（(,)kjkjuuxttt）：Otutxutxukjkjkj,,12112,,OtutxutxukjkjkjhOxuhtxutxukjkjkj,,1hOxuhtxutxukjkjkj,,12112,,hOxuhtxutxukjkjkj222211,,2,hOxuhtxutxutxukjkjkjkj可得到以下几种最简差分格式（一）向前差分格式14.1kjkjuu1jkjkjkjfhuuua2112jjxff24.1jjjxu0，ku0=kNu=0其中1,,1,0Nj，1,,1,0Mk。取2har为网比，则进一步有14.11kju=kjru1+r21kju+kjru1+jf此差分格式是按层计算：首先，令0k，得到1ju=01jru+r210ju+01jru+jf于是，利用初值jjjxu0和边值ku0=kNu=0，可算出第一层的1ju，1,,1,0Nj。再由14.1取1k，可利用1ju和ku0=kNu=0算出2ju，1,,1,0Nj。如此下去，即可逐层算出所有kju（1,,1,0Nj，1,,1,0Mk）。由于第1k层值可以通过第k层值直接得到，如此的格式称为显格式。并视kju为kjtxu,的近似值。若记TkNkkkuuu121,,,u，TNxxx121,,,，TNxfxfxf121,,,f则显格式14.1可写成向量形式011,,1,0,ufAuuMkkk其中rrrrrrrrrr21002100210021A若记22xuatuLukjkjkjhuuuL112112huuuakjkjkj那末截断误差（1.5）uRkjkjkjhLutxuL,1=Otxturkj)~,~(2112122=2hO其中(,)jkxt是矩形11jjxxx，1jkttt中某一点。事实上，uRkjkjxu222+2Okjxuha442ˆ12=kjxu222+2O22222ˆ112Otuahakj=211212ah222~Otukj=21121r222~Otukj=2hO。这里222244xuxaxuatuaxa122tux22txu2322tutut222xuattxu23故22tu44244xuaxuaa，从而44xu221tua（二）向后差分格式16.1kjkjuu1jkjkjkjfhuuua2111112jjxff26.1jjjxu0，ku0=kNu=0其中1,,1,0Nj，1,,1,0Mk。取2har为网比，则进一步有16.1rkju1+r211kjur11kju=kju+jf按层计算：首先，取0k，则利用初值jjjxu0和边值ku0=kNu=0，来确定出第一层的1ju，1,,1,0Nj，即求解方程组：r11ju+r211jur11ju=0ju+jf1,,1,0Nj，ku0=kNu=0。求出1ju，在由14.1取1k，可利用1ju，解出2ju，1,,1,0Nj。如此下去，即可逐层算出所有kju，1,,1,0Mk。如此每层必须解一个三对角线性方程组的格式称为隐格式。并视kju为kjtxu,的近似值。直观地说，采用显式格式进行求解既方便又省工作量。但是，后面我们将看到，有些情况用隐式格式更为便利。1.2.3Grank-Nicholson法将向前差分格式和向后差分格式做算术平均，得到的差分格式称之为六点对称格式，也称为Crank-Nicholson格式：18.1kjkjuu1jkjkjkjkjkjkjfhuuuhuuua211111211222jjxff28.1jjjxu0，ku0=kNu=0进一步，18.12r11kju+r11kju2r11kju=2rkju1+r1kju2rkju1+jf按层计算：首先，取0k，则利用初值jjjxu0和边值ku0=kNu=0，来确定出第一层的1ju，1,,1,0Nj，即求解方程组：2r11ju+r11ju2r11ju=2r01ju+r10ju2r01ju+jf1,,1,0Nj，ku0=kNu=0。求出1ju，在由18.1,取1k，可利用1ju，解出2ju，1,,1,0Nj。如此下去，即可逐层算出所有kju，1,,1,0Mk。若记22xuatuLukjkjkjhuuuL13jkjkjkjkjkjkjfhuuuhuuua211111211222在1(,)(,()/2)jkkxtxtt处作Taylor展开，可以算出截断误差为（1.7）uRkjkjkjhLutxuL,3=22hO。+（四）Richardson格式（1.10）211kjkjuu2112huuuakjkjkj+jf进一步110.11kju=r2（kju1kju2+kju1）+1kju+2jf这是三层显式差分格式。显然截断误差的阶为22hO。为使计算能够逐层进行，除初值0ju外，还要用到1ju。它可以用其他双层格式提供。Richardson格式的矩阵形式为：另算10111,,1,uufuuCuMkkkk其中21001210012100122rC2稳定性与收敛性抛物方程的两层差分格式可以统一写成向量形式:(2.1)1kkAUBUF其中1111(,),(,,)kkkTNNUuuFff，A和B是1N阶矩阵。我们假定A可逆，即（2.1）是唯一可解的。对于显格式，A等于单位矩阵I。三层格式可以通过引入新变量1kkkUWU化成两层格式。假设差分解的初始值（其实可以是任一层的值）有误差，以后各层计算没有误差，让我们来考察初始误差对以后各层的影响。令kU和kV分别是以0U和0V为初始值由差分格式（2.1）得到的两组差分解，则kkkWUV满足（2.2）1kkAWBW因此，按初值稳定应该意味着0kWKW。这就导致如下定义：假设0F，我们称差分格式（2.1）按初值稳定，如果存在正常数0和K，使得以下不等式成立：（2.2）0kUKU，010,0,0NURkT这里是1NR上的某一个范数，例如1/2121NjjUhu类似地，假设00U，我们称差分格式（2.1）按右端稳定，如果存在正常数0和K，使得以下不等式成立：（2.2）kUKF，00,0kT可以证明，差分格式若按初值稳定，则一定按右端稳定。因此，这时我们简单地称差分格式稳定。前面讨论的向前差分格式（1.4）当网比12r时稳定，当12r时不稳定。这就意味着给定空间步长h以后，时间步长必须足够小，才能保证稳定。而向后差分格式（1.6）和Grank-Nicholson格式（1.8）则对任何网比r都是稳定的，时间步长可以取得大一些，从而提高运算效率。Richardson格式则对任意网比都是不稳定的。因此，虽然Richardson格式是个显格式，截断误差又很小，但是却不可用。如果某个差分格式的截断误差当和h趋于0时随之趋于0，则称这个差分格式是相容的。可以证明：若差分格式是相容的和稳定的，则它是收敛的，并且差分解与微分解之间误差的阶等于截断误差的阶。因此，当网比12r时，向前差分格式（1.4）有收敛阶2()Oh。对任何网比，向后差分格式（1.6）有收敛阶2()Oh，而Crank-Nicholson格式（1.8）有收敛阶22()Oh。3．高维抛物方程差分法考虑如下二维抛物方程的差分格式。（3.1）2222,,(0,),0(,,0)(,)(0,,)(,,)(,0,)(,,)0uuuxylttxyuxyxyuytulytuxtuxlt取空间步长/hlN，时间步长0。作两族平行与坐标轴的网线jxxjh，kyykh，其中,0,1,,jkN，将矩形区域(0,)(0,)ll分割成2N个小矩形。记njku为网格节点(,,)jknxyt上的差分解。前述各种一维差分格式都可以直接用于以（3.1）为代表的二维以至更高维的抛物方程。例如，向前差分格式成为（3.2）11,1,,1,12222nnnnnnnnjkjkjkjkjkjkjkjkuuuuuuuuhh0,1,,1;,1,,1TnjkN实际计算时，先令0n，利用已知的0(,)jkjkuxy等等，对,1,,1jkN，用（3.2）算出1jku。而由边值条件，补充得到11110,,,0,0kNkjjNuuuu。下一步，令1n，利用已知的第1层的差