有限差分法求解抛物型方程说明

有限差分法求解抛物型方程偏微分方程只是在一些特殊情况下，才能求得定解问题解的解析式，对比较复杂的问题要找到解的解析表达式是困难的，因此需采用数值方法来求解.有限差分法是一种发展较早且比较成熟的数值求解方法，只适用于几何形状规则的结构化网格.它在微分方程中用差商代替偏导数，得到相应的差分方程，通过解差分方程得到微分方程解的近似值.本章主要介绍有限差分法的基本思想，并给出一些具体的数值实例.§1差分方法的基本思想有限差分法把偏微分方程的求解区域划分为有限个网格节点组成的网格，主要采用Taylor级数展开等方法，在每个网格节点上用有限差分近似公式代替方程中的导数，从而建立以网格节点上的函数值为未知数的代数方程组.有限差分格式，从格式的精度来划分，有一阶格式、二阶格式和高阶格式.从差分的空间形式来考虑，可分为中心格式和逆风格式.考虑时间因子的影响，差分格式还可以分为显格式、隐格式和显隐交替格式等.目前常见的差分格式，主要是上述几种格式的组合，不同的组合构成不同的差分格式.泰勒级数展开法对有限差分格式的分类和公式的建立起着十分重要的作用.下面采用泰勒展开式导出一个自变量系统的若干有限差分表达式.首先考虑单变量函数()ux，如图1把区域x离散为一批结点，记0()(),=0,1,2,iiuxuxihui图1单变量函数离散化函数()ux在点ix处的泰勒展开式为23()()()()()2!3!iiiiiuxuxuxhuxuxhhh（1）或23()()()()()2!3!iiiiiuxuxuxhuxuxhhh（2）式(1)和(2)重新整理可得2()()()()()2!3!iiiiiuxhuxuxuxuxhhh（3）和2()()()()()2!3!iiiiiuxuxhuxuxuxhhh（4）于是给出在点ix处函数u的一阶导数的两个近似公式1()()()iiiiiuxhuxuuuxhh（5）1()()()iiiiiuxuxhuuuxhh（6）因为级数被截断，这两个近似公式肯定要产生误差，此误差与h同阶，形式分别为()(),,2()(),.2iiiiiihEuOhxxhhEuOhxhx若把式(3)和(4)相加并求()iux，可得11()()()22iiiiiuxhuxhuuuxhh（7）其截断误差与2h同阶，形式为22()(),,6iiihEuOhxhxh若把式(3)和(4)相减并求()iux，可得1122()2()()2()iiiiiiiuxhuxuxhuuuuxhh（8）其截断误差与2h同阶，其形式为22()(),,12iiihEuOhxhxh我们可继续用这种方式来推导更复杂的公式，类似的公式还有很多，这里不再一一列举.公式（5）、（6）分别称为一阶向前、向后差分格式，这两种格式具有一阶计算精度，公式（7）、（8）分别称为一阶、二阶中心差分格式，这两种格式具有二阶计算精度.图2二维区域网格剖分上面的结果可直接推广使用于导出二元函数(,)uxy的许多有限差分近似公式.如图7.2，把求解区域进行网格剖分，使12(,)(,),,=0,1,2,ijijuxyuihjhuij其中x方向的网格间距为1,hy方向的网格间距为2,h整数i和j分别表示函数(,)uxy沿x坐标和y坐标的位置.二元函数(,)uxy对x求偏导时y保持不变，对y求偏导时x保持不变，根据向前差分公式（7.5）可以给出在点(,)ijxy处函数(,)uxy的一阶偏导数的两个近似公式1,,1(,)ijijijuxyuuxh（9）,1,2(,)ijijijuxyuuyh（10）相类似地，根据二阶中心差分格式（8）可以得到函数(,)uxy的二阶偏导数的近似公式21,,1,221(,)2ijijijijuxyuuuxh（11）2,1,,1222(,)2ijijijijuxyuuuyh（12）下面我们推导函数(,)uxy的二阶混合偏导数2uxy在(,)ijxy的有限差分表达式.根据一阶中心差分格式（7），112111,11,11,11,122121221,11,1(,)(,)(,)1()21()()222ijijijijijijijijijiuxyuxyuxyOhxyhyyuuuuOhOhhhhuuu1,11,1124jijuhh二维有限差分近似可以直接推广到三维空间或三维空间加一维时间的情形.定义1当步长趋于零时，差分方程的截断误差趋于零，则称差分格式与微分方程是相容的.定义2当步长趋于零时，差分方程的解收敛于微分方程的解，则称差分格式是收敛的.定义3当差分方程的解由于舍入误差的影响，所产生的偏差可以得到控制时，则称差分格式是稳定的.§2抛物型方程的有限的差分法为了说明如何使用有限差分法来求解偏微分方程，本节我们给出以下几个数值实例.算例1考虑一维非齐次热传导方程的初边值问题：2212(,),01,01,(,0)(),01,(0,)(),(1,)(),01.uuafxtxttxuxqxxutgtutgtt(7.13),其中2,a函数11(,)[cos()2sin()],22xfxtett初始条件1()sin,2xqxe左、右边界条件分别为11()sin(),2gtt21()sin()2gtet.该定解问题的解析解为1(,)sin(),(,)[0,1][0,1].2xuxtetxt将求解区域{(,)|,0}xtaxbtT进行网格剖分，[,]ab作m等分，[0,]T作n等分，记,,baThmn则,0,,0ikxaihiMtkkn对该问题建立如下向前差分格式：11122,11,11,kkkkkkiiiiiiuuuuuafimknh（14）(,0)(),1,iiuxqxim（15）12(,)(),(,)(),1.kkkkuatgtubtgtkn（16）令2rah，差分格式（7.14）整理得111(12),11,11.kkkkkiiiiiurururufimkn（17）显然时间在1kt上的每个逼近值可独立地由kt层上的值求出。该格式为显格式，采用显格式时，应注意时间步长和空间步长的选取，当12r时向前差分格式是稳定的.我们采用步长0.1h和0.0025，选取0.25t时的数据进行比较，得到精确解与近似解的最大误差是9.00670e-005.表1算例10.25t时节点处数值解、精确解和误差的绝对值（显式格式）节点0.00.10.20.30.4近似解0.2474040.2734520.3022320.3340310.369167精确解0.2474040.2734240.3021800.3339600.369083绝对误差（10-4）00.2835440.5229660.7104910.838755节点0.50.60.70.80.9近似解0.4079900.4508880.4982900.5506700.608551精确解0.4079000.4507990.4982100.5506080.608516绝对误差（10-4）0.9006700.8896480.7994350.6243880.359370数值解与精确解的比较见下图：图3t=0.25精确解与近似解的比较图4t=0.25精确解与近似解的绝对误差对该问题还可建立如下向后差分格式：11122,11,11,kkkkkkiiiiiiuuuuuafimknh（18）差分格式（14）整理得111(12),11,11.kkkkkiiiiirururuufimkn（19）显然时间在kt层上的逼近值需要通过求解一个三对角线性方程组得到。该格式为隐格式，该隐格式对于任意网格比r都是稳定的.我们采用步长0.1h和0.0025，选取0.25t时的数据进行比较，得到精确解与近似解的最大误差是1.274828e-005.数值解与精确解的比较结果如下图：图5t=0.25精确解与近似解的比较图6t=0.25精确解与近似解的绝对误差算例2考虑一维非线性Chaffee-Infante方程的初边值问题：11(())(),(,)[0,1][0,1],(,0)(),[0,1],(,)(),(,)(),[0,1].uukxfuxttxxuxqxxuatgtubtgtt（20）其中扩散系数()1,kx3(),fuuu初始条件1111()tanh(1),2222qxx左右边界条件分别为：111113()tanh(1),22224gtat211113()tanh(1)22224gtbt.该问题的解析解为11113(,)tanh(1),(,)[,][0,].22224uxtxtxtabT.将求解区域进行剖分，方法同上，对该问题建立如下向前差分格式：1111122(),11,11,kkkkkkkkiiiiiiiiiiiuukkuuuuukfuimknhhh(,0)(),1,iiuxqxim12(,)(),(,)(),1.kkkkuatgtubtgtkn采用步长0.1h和0.005,我们选取0.5t时的数据进行比较，得到精确解与近似解的最大误差是6.866286e-005.表2算例20.5t时节点处数值解、精确解和误差的绝对值节点0.00.10.20.30.4近似解0.2227000.2106750.1991360.1880830.177512精确解0.2227000.2107000.1991810.1881410.177578绝对误差（10-4）00.2517250.4463340.5834810.663365节点0.50.60.70.80.9近似解0.1674170.1577910.1486250.1399100.131632精确解0.1674860.1578560.1486820.1399520.131656绝对误差（10-4）0.6866280.6543120.5677560.4285880.238650运行结果如下图，图7t=0.5精确解与近似解的比较图8t=0.5精确解与近似解的绝对误差