§4 偏微分方程的数值解法

一、 差分法

差分法是常用的一种数值解法.它是在微分方程中用差商代替偏导数，得到相应的差分方程，通过解差分方程得到微分方程解的近似值. 1. 网格与差商

在平面 (x,y)上的一以S为边界的有界区域D上考虑定解问题.为了用差分法求解，分别作平行于x轴和y轴的直线族.

?x?xi?ih (i,j=0,?1,?2,…,?n) ?y?y?jhi?作成一个正方形网格，这里h为事先指定的正数，称为步

长；网格的交点称为节点，简记为(i,j).取一些与边界S接近的网格节点，用它们连成折线Sh，Sh所围成的区域记作Dh.称Dh内的节点为内节点，位于Sh上的节点称为边界节点（图14.7）.下面都在网格Dh + Sh上考虑问题：寻求各个节点上解的近似值.在边界节点上取与它最接近的边界点上的边值作为解的近似值，而在内节点上，用以下的差商代替偏导数：

?u1??u?x?h,y??u?x,y???xh?u1??u?x,y?h??u?x,y???yh

图14.7

?2u1?u?x?h,y??2u?x,y??u?x?h,y????x2h2?2u1?u?x,y?h??2u?x,y??u?x,y?h????y2h2

?2u1?2?u?x,y?h??u?x?h,y??u(x,y?h)?u?x,y???x?yh11 注意， 1? 式中的差商?u?x?h,y??u?x,y??称为向后差商，而?u?x,y??u?x?h,y??称为向

hh1前差商，?u?x?h,y??u?x?h,y??称为中心差商.也可用向前差商或中心差商代替一阶偏导数.

2h 2? x轴与y轴也可分别采用不同的步长h,l，即用直线族

?x?xi?ih (i,j=0, ±1, ±2,? ) ?y?y?jhj?作一个矩形网格.

2. 椭圆型方程的差分方法

[五点格式] 考虑拉普拉斯方程的第一边值问题

??2u?2u?x,y?D???0???x2?y2 ??u???x,y?S??式中?(x,y)为定义在D的边界S上的已知函数.

采用正方形网格，记u(xi,yj)=uij ，在节点(i,j)上分别用差商

ui?1,j?2uij?ui?1,jui,j?1?2uij?ui,j?1 ,h2h2?2u?2u代替2,2，对应的差分方程为

?x?yui?1,j?2uij?ui?1,jui,j?1?2uij?ui,j?1??0 （1） 22hh或

1uij?ui?1,j?ui?1,j?ui,j?1?ui,j?1

4即任一节点(i,j)上uij的值等于周围相邻节点上解的值的算术平均，这种形式的差分方程称为五点格式，在边界节点上取

uij??xi*,y*j??i,j??Sh? （2）

??式中(xi*,yj*)是与节点(i,j)最接近的S上的点.于是得到了以所有内节点上的uij值为未知量的若干个

线性代数方程，由于每一个节点都可列出一个方程，所以未知量的个数与方程的个数都等于节点的总数，于是，可用通常的方法（如高斯消去法）解此线性代数方程组，但当步长不很大时，用高斯消去法将会遇到很大困难，可用下面介绍的其他方法求解.

若h?0时，差分方程的解收敛于微分方程的解，则称差分方程为收敛的.

在计算过程中，由于进行四则运算引起舍入误差，每一步计算的舍入误差都会影响以后的计算结果，如果这种影响所产生的计算偏差可以控制，而不至于随着计算次数的增加而无限增大，则称差分方程是稳定的.

[迭代法解差分方程] 在五点格式的差分方程中，任意取一组初值{uij}，只要求它们在边界节点(i,j)上取以已知值?(xi*,yj*)，然后用逐次逼近法（也称迭代法）解五点格式：

1n??n?1??n??n??n??n?0,1,2,?? uij?ui??1,j?ui?1,j?ui,j?1?ui,j?14逐次求出{uij(n)}.当(i+1,j)，(i－1,j)，(i,j－1)，(i,j+1)中有一点是边界节点时，每次迭代时，都要在这一点上取最接近的边界点的值.当n?∞时，uij(n)收敛于差分方程的解，因此n充分大时，{uij(n)}可作差分方程的近似解，迭代次数越多，近似解越接近差分方程的解.

[用调节余数法求节点上解的近似值] 以差商代替Δu时，用节点(i+1,j)，(i－1,j)，(i,j+1)，(i,j－1)上u的近似值来表示u在节点(i,j)的值将产生的误差，称此误差为余数Rij，即

u?xi?h,yj??u?xi?h,yj??u?xi,yj?h??u?xi,yj?h??4u?xi,yj??Rij

?? 设在(i,j)上给uij以改变量?uij，从上式可见Rij将减少4?uij，而其余含有u(xi,yj)的差分方程中的余数将增加?uij，多次调整?uij的值就可将余数调整到许可的有效数字的范围内，这样可获得各节点上u(x,y)的近似值.这种方法比较简单，特别在对称区域中计算更简捷.

例 求Δu=0在内节点A,B,C,D上解的近似值.设在边界节点1,2,3,4上分别取值为1,2,3,4（图14.8）

解 记u(A)=uA，点A,B,C,D的余数分别为

图14.8

－4uA+ uB+ uc +5=RA uA－4 uB + uD+7=RB

u －4 uc+ uD+3=RC A

uB+ uc－4uD+5=RD

以边界节点的边值的算术平均值作为初次近似值，即

uA(0)=uB(0)=uC(0)=uD(0)=2.5

则相应的余数为：

RA=0, RB=2, RC= －2, RD=0

最大余数为±2.先用δuC=－0.5把RC缩减为零，uC相应地变为2，这时RA, RD也同时缩减(－0.5)，新余数是RA=－0.5,RB=2,RC?0, RD=－0.5.类似地再变更δuB=0.5，从而 uB变为3，则得新余数为RA?RB?RC?RD?0.这样便可消去各节点的余数，于是u在各节点的近似值为：

uA=2.5, uB=3, uC=2, uD=2.5

现将各次近似值及余数列表如下：

次调 整 值 数 0 1 δuC = －0.5 2 δuB = 0.5 结果近似值

uA 2.5 2.5 2.5 2.5 RA 0 －0.5 0 uB 2.5 2.5 3 3 第n次近似值及余数 RB uC RC 2 2.5 －2 2 2 0 0 2 0 2 uD 2.5 2.5 2.5 2.5 RD 0 －0.5 0 [解重调和方程的差分方法] 在矩形D(x0≤x≤x0+a,y0≤y≤y0+a)中考虑重调和方程

?4u?4u?4u2?u?4?222?4?0

?x?x?y?ya取步长h?，引直线族

n?x?x0?ih (i, j = 0, 1, 2,?, n) ??y?y0?jh作成一个正方形网格.用差商代替偏导数

1u?x,y???8?u?x?h,y??u?x?h,y??u?x,y?h??u?x,y?h??20?2?u?x?h,y?h??u?x?h,y?h??u?x?h,y?h??u?x?h,y?h?? ??u?x?2h,y??u?x?2h,y??u?x,y?2h??u?x,y?2h???上式表明了以(x,y)为中心时，u(x,y)的函数值与周围各点函数值的关系，但对于邻近边界节点的点(x,y)，如图14.9中的A，就不能直接使用上式，此时将划分网格的直线族延伸，在延伸线上定出与边界距离为h的点，称这些点为外邻边界节点，如图14.9以A为中心时，点E,C为边界节点，点J,K为E,C的外邻边界节点，用下法补充定义外邻边界节点J处函数的近似值uJ，便可应用上面的公式.

1? 边界条件为

?uuS??1?P?,??2?P??P?S? 图14.9

?xS