第10章 偏微分方程数值解_2

10.3抛物型方程的差分解法抛物型方程是指如下形式的方程：),,(222222zyxfzuyuxuatu很多实际的物理问题都可以用这类方程描述：热传导方程：),(22txfxuatu现以热传导方程为例，介绍抛物型方程的有限差分格式。设热传导方程：定解条件),(22txfxuatu(1)0),1(),0()(0tutuxut(2)求（1）满足（2）的解。10.3.1矩形网格用两组平行直线族xj=jh,tk=k(j=0,1,…，k=0，1,…)构成的矩形网覆盖了xt平面，网格点（xj,tk）称为结点，简记为（j,k），h、为常数，分别称为空间步长及时间步长，或称h为沿x方向的步长，称为沿t方向的步长，，N为正整数。在t=0上的结点称为边界结点，其余所有属于内的结点称为内部结点。Ttxc0,txoh(xj,tk)10.3.2．古典差分格式于平面区域上考虑传导方程：TtxtxG0,10),(),(22txfxuatuLUa为正常数（3）Tttutuxxxu00),1(),0(10)()0,((4)于结点（j,k）处偏导数与差商之间有如下近似的关系：)(),(),(1otutxutxukjkjkj)(),(),(2)(222211hoxuhtxutxutxukjkjkjkj利用上述表达式得到LU在(j,k)处的关系式：2111),(),(2)(),(),(htxutxutxuatxutxukjkjkjkjkj)(][2hoLUkj(5)),(kjkjtxff视为u(xj,tk)的近似值。kju令，j=1,2,…,N–1；k=0,1,2,…kjkjkjkjkjkjfhuuuauu21112则有：（6）差分方程（6）称为解热传导方程（3）的古典显格式，它所用到的结点如下图：****（j,k）将（6）写成便于计算的格式：kjkjkjkjkjkjfuuuruu1112（7）2/har称为网比，利用（7）及初边值条件（4）在网格上的值0)(00kNkjjjuuxu（8）即可算出k=1,2,…，各层上的值。kju截断误差阶为0(+h2)。为了提高截断误差的阶，可以利用中心差商：)(2),(),(211otutxutxukjkjkjj=1,2,…,N–1；k=0,1,2,…kjkjkjkjkjkjfhuuuauu2111122（9）得到Richardson格式，其结点图为：****（j,k）*截断误差阶为o(2+h2)，较古典显格式高。将（9）式改写成适于计算的形式：j=1,2,…,N–1；k=1,2,…r=a/h2称为网比，（10）式中出现了三层网格上的值，kjkjkjkjkjkjfuuuuru2221111（10）才能逐层计算。故需要事先求得第k－1层的值1kjukju和第k层的值,如果利用向后差商)(),(),(1oxutxutxukjkjkjj=1,2,…,N–1；k=0,1,2,…kjkjkjkjkjkjfhuuuauu21112(11)kjkjkjkjkjkjfuuuuru1112(12)j=1,2,…,N–1；k=0,1,2,…古典隐格式，其结点图为：（j,k）****截断误差为o(+h2)，与古典显格式相同。10.3.3．六点对称格式取该点的中心差商，从而对于方程（3）式，在点列方程，，21,kjtu22xu)(212/1ouuxukjkjkj)(22/12/12/12/12/122hohxuxuxukjkjkj)(122/112/12/12/11hohuuhuuhkjkjkjkjkjkjkjuuu1112/1121kjkjkjuuu12/121kjkjkjfff12/121将以上各式代入（3）式得到差分方程：2/12/112/12/11212kjkjkjkjkjkjfuuuhauu2/1111111122222kjkjkjkjkjkjkjfuuuuuuha2/1111111121222kjkjkjkjkjkjkjkjkjfuuuuuuhauu111112)1(2kjkjkjururur2/1112)1(2kjkjkjkjfururur整理,得此即六点对称格式，也称为Crank-Nicolson格式，所用结点图为：***k+1***kj+1jj–1（13）10.3.4．稳定性（1）当步长无限缩小时，差分方程的解是否逼近于微分方程（2）计算过程中产生的误差在以后的计算中是无限增加，还是可以控制？（稳定性）的解？（收敛性）稳定性问题是研究抛物型差分方程的一个中心课题！考察Richardson格式的稳定性。用表示计算所产生的误差，如果右端无误差存在，kje1kjukjf则满足：kje111122kjkjkjkjkjeeeere取21r11112kjkjkjkjkjeeeee（14）假设k-1层之前无误差存在。即，而在第k层产生了01kje误差。，这一层其它点也无误差，而且在计算过程kje0中不再产生新的误差，利用（14）式算出误差的传播如下表：r=½时Richardson格式的误差传播jj0–4j0–3j0–2j0–1j0j0+1j0+2j0+3j0+4k-2-474-617-2417-6-831-6889-6831-8-1049-144277-388277-14449-1071-260641-1091311-109641-26071r≤1/2时古典显格式的误差传播jj0–4j0–3j0–2j0–1j0j0+1j0+2j0+3j0+4k0.500.50.2500.500.250.12500.37500.37500.1250.062500.2500.37500.2500.0625如果选用r=½时的古典显格式，误差方程为：kjkjkjeee11121差分格式关于初值稳定的实际含义是：如果其解在某一层存在误差，则由它引起的以后各层上的误差不超过原始误差的M倍（M为与无关的常数）。因此，在稳定的条件下，只要初始误差足够小，以后各层的误差也能足够小。以上构造的几种差分格式中，古典显格式：r≤1/2时稳定古典隐格式：绝对稳定Richardson格式：绝对不稳定六点对称格式：绝对稳定。稳定性概念：（1）水流为稳态和无其反应情况下的溶质运移方程xCvxCDtC22考虑六点对称格式，在点列方程，21,kj)(212/1oCCtCkjkjkjthx令，，有：2/112/12/1122/1222kjkjkjkjCCChDxC)(22/12/12/12/12/1hohCCvxCvkjkjkj)(2222211111112hoCCCCCChDkjkjkjkjkjkj)(2222/112/12/12/11hohCCCCvkjkjkjkj)(222/112/11hoCChvkjkj)(22111111hoCCCChvkjkjkjkj舍去o(h2)，o(2)得到六点对称格式：kjkjkjkjkjkjkjkjCCCCCChDCC11111112122221122hDDhvv令上式变为1111111111)()21()(kjkjkjCDvCDCDvkjkjkjkjCCCChv1111112kjkjkjCvDCDCDv1111111)()21()(kjkjkjjCvDCDCvDH1111111)()21()()(),21(),(11111DvFDBDvAjjj令上式可以写成：jjjjjjjHCFCBCA11j=1,2,…,N–1对时间变量用向后差分，对空间变量用中心差分，hCCvhCCCDCCkjkjkjkjkjkjkj22112111可得到隐格式：hvvhDD121,令整理得：1111111)()21()(kjkjkjkjCCvDCrCvD（2）．非稳态方程xqCzCvDztCsh)(),()(非稳态一维垂直流情况下，土壤溶质的基本方程为：式中容积含水量的求法如下：zhKzhhKzthhC)()()(式中，h为负压水头；t为时间；z为到原点的距离（cm），先解非饱和垂直水流方程dhdhC)(向下为正；C(h)为容水度，，K(h)为土壤导水率，可由一些常用的经验公式算出。求出h后，用水分特征曲线换算成相应的值。v=q/，q为通量。)(2111122/12/12/12/111kikikikikikikikikikiCCCCzGDtCC8)(212/12/12/12/12/12/1kikikikikikitvvqzG式中：2/12/111112/12/12/12/12kikikikikikikikzDq)(2111122/12/12/12/1kikikikikikiCCCCzGDzCCqCCqkikikikikiki2)()(1112/12/112/12/1)(2111112/12/1kikikikiki2212/112/1kikikikikikiikiikiikiiHCFCBCA11111将以上方程整理后可写成：2/12/12/12/12/12/1kikikiizqDGA2/12/12/12/12/12/12/12/12/1/kikikikikizqDGDrB2/12/12/12/1kikiiDGFkikikikikikikikiikiiiCrGDGDqzCFCAH2/12/12/12/12/12/12/12/12/12/11122ztr2、有限元法设有微分方程定义在由边界围成的区域以上，L为微分算子。设{j}（j=1,2,…，n,…），是一完备的函数系。伽辽金方法是求形如的近似解，其中aj(j=1,2,…,n)为待定常数。un称为试探函数，j称为形状函数（或基函数，插值函数，为{j}(j=1,2,…,n,…)中前n个线性无关的函数）。fLu（1）0fLujnjjnau1(2)若u是方程（1）的精确解，则必有在Lu和f是连续函数的条件下，就等价于,,,2,10)(nidfLui