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

 

一、    差分法

 

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

    1.  网格与差商

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

      (i,j=0,1,2,…,n) 

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

    注意， 1°  式中的差商称为向后差商，而称为向前差商，称为中心差商.也可用向前差商或中心差商代替一阶偏导数.

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

   (i,j=0, ±1, ±2 )

作一个矩形网格.

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

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

式中(x,y)为定义在D的边界S上的已知函数.

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

代替，对应的差分方程为

                                     （1）

或

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

                                                           （2）

式中(x_i*,y_j*)是与节点(i,j)最接近的S上的点.于是得到了以所有内节点上的u_ij值为未知量的若干个线性代数方程，由于每一个节点都可列出一个方程，所以未知量的个数与方程的个数都等于节点的总数，于是，可用通常的方法（如高斯消去法）解此线性代数方程组，但当步长不很大时，用高斯消去法将会遇到很大困难，可用下面介绍的其他方法求解.

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

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

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

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

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

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

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

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

－4u_A+  u_B+  u_c        +5=R_A

   u_A－4 u_B         + u_D+7=R_B

       u_A         －4 u_c+ u_D+3=R_C

             u_B+  u_c－4u_D+5=R_D

 

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

u_A⁽⁰⁾=u_B⁽⁰⁾=u_C⁽⁰⁾=u_D⁽⁰⁾=2.5

则相应的余数为：

R_A=0,  R_B=2,  R_C= －2,  R_D=0

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

u_A=2.5,  u_B=3,  u_C=2,  u_D=2.5

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

 

次数	调 整 值	第n次近似值及余数
		uA	RA	uB	RB	uC	RC	uD	RD
0 1 2	δuC = －0.5 δuB =  0.5	2.5 2.5 2.5	0 －0.5 0	2.5 2.5 3	2      2 0	2.5 2 2	－2 0 0	2.5 2.5 2.5	0 －0.5 0
结果近似值		2.5		3		2		2.5

 

    [解重调和方程的差分方法]  在矩形D(x₀≤x≤x₀+a,y₀≤y≤y₀+a)中考虑重调和方程

取步长，引直线族

    (i, j = 0, 1, 2 n)

作成一个正方形网格.用差商代替偏导数

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

1°  边界条件为

时，定义u_J=u_A－2₂(E)h.

2°    边界条件为

时，定义u_J=2₁(E)－u_A－h²₂(E).

    [其他与Δu有关的网格]

    1°  三角网格（图14.10(a)）

    取P₀(x,y)为中心，它的周围6个邻近节点分别为：

则                                             

式中u_i=u(P_i), u₀=u(P₀)，R表示余项.

    2°  六角网格（图14.10(b)）
    取P₀(x,y)为中心，它的三个邻近节点分别为

则                           .

图14.10

 

    3°  极坐标系中的网格（图14.10(c)）

    取P₀(r,)为中心，它的四个邻近节点分别为

而拉普拉斯方程

的相应的差分方程为

    3.  抛物型方程的差分方法

    考虑热传导方程的边值问题

将[0,b]分为n等份，每段长为.引两族平行线（图14.11）

图14.11

 

x=x_i=ix      (i=0,1,2n)

y=y_j=jt      (j=0,1,2  t 取值见后)

作成一个长方形的网格，记u(x_i,t_j)为u_ij，节点(x_i,t_j)为(i,j)，在节点(i,j)上分别用

代替，于是边值问题化为差分方程

记，差分方程可写成

                      （1）

由此可按t增加的方向逐排求解.在第0排上u_i0的值由初值(ix)确定，j+1排u_i,j+1的值可由第j排的三点(i+1,j),(i,j),(i－1,j)上的值u_i+1,j, u_ij,u_i-1,j确定，而u_0,j+1,u_n,j+1已由边界条件₁((j+1)t)及₂((j+1)t)给定，于是可逐排计算一切节点上的u_ij值.当(x), ₁(x)和₂(x)充分光滑，且时，差分方程收敛而且稳定.所以利用差分方程（1）计算时，必须使，即.

    热传导方程还可用差分方程

代替，此时如已知前j排u_ij的值，为求第j+1排的u_i,j+1 必须解包含n－1个未知量的线性代数方程组，这种差分方程称为隐式格式的差分方程，前面所提的差分方程称为显式格式差分方程.

    隐式格式差分方程对任意的λ都是稳定的.

    4.  双曲型方程的差分方法

    考虑弦振动方程的第一边值问题

    用矩形网格，列出对应的差分方程：

    记与上段一样，利用和在第0排及第1排的已知数值（初始条件）u_i0 , u_i1可计算u_i2，然后用已知的u_i1 , u_i2及可计算u_i3,类似地可确定一切节点上的u_ij值.

    当(x),(x),₁(x)和₂(x)充分光滑，且ω≤1时，差分方程收敛且稳定，所以要取.

 

二、    变分方法

 

    1.  自共轭边值问题

    将§3定义的共轭微分算子的概念推广到一般方程.

    设D是中的有界区域，S为其边界，在上考虑2k阶线性微分方程

的齐次边值问题

式中f(x)是D内的已知函数，l_ju是线性微分算子.

    将 分部积分k次得

式中(u,v)是一个D上的积分，其被积函数包含u,v的k阶导数；R_j和是定义在边界S上的两个线性微分算子.再将(u,v)分部积分k次得

式中L*是一个2k阶的微分算子，称为L的共轭微分算子.若L=L*，则称L为自共轭微分算子.从上面可推出格林公式

如从l_ju|_S=l_jv|_S=0可推出在边界S上

则称l_ju|_S=0为自共轭边界条件.如果微分算子及边界条件都是自共轭的，则称相应的边值问题为自共轭边值问题，此时有

    每个边值问题对应于某希尔伯特空间H（例如L₂(D)，见第九章§7）中的一个算子A，其定义域M_A 是H中一线性稠密集合，它由足够次连续可微且满足边界条件的函数组成，在M_A上，Au的数值与Lu的数值相同，从而求解边值问题化为解算子方程

的问题.

    设A为定义在实的希尔伯特空间H中的某线性稠密集合M_A上的线性算子.若对于M_A的任意非零元素成立

(Au,v)=(u,Av)

则称A为对称算子.若对任意非零元素u成立

则称A为正算子.如成立更强的不等式

(Au,u)≥r||u||²      (r>0)

则称A为正定算子.此处(u,v)表示希尔伯特空间的内积，||u||²=(u,u).

    2.  变分原理与广义解

    定理  设A是正定算子，u是方程Au=f在M_A上的解的充分必要条件是： u使泛函

F(u)=(Au,u)－2(f,u)

取极小值.

    上述将边值问题化为等价的求泛函极值问题的方法称为能量法.在算子的定义域不够大时，泛函F(u)的极值问题可能无解.不过对于正定算子，可以开拓集合M_A，使在开拓了的集合上，泛函的极值问题有解.为开拓M_A，在M_A上引进新的内积[u,v]=(Au,v)，定义模||u||²=[u,u]=(Au,u)，在模||u||的意义下，补充极限元素，得到一个新的完备希尔伯特空间H₀，在H₀上，泛函F(u)仍然有意义，而泛函的极值问题有解.但必须注意，此时使泛函F(u)取极小的元素u₀不一定属于M_A，因此它不一定在原来的意义下满足方程Au=f及边界条件.称u₀为广义解.

    3.  极小化序列与里兹方法

    在处理变分问题中，极小化序列起着重要的作用.考虑泛函

F(u)=(Au,u)－2(f,u)

以d表示泛函的极小值.设在希尔伯特空间中存在一列元素{u_n}  (n=1,2)，使

则称{u_n}为极小化序列.

    定理  若算子A是正定的，则F(u)的每一个极小化序列既按H空间的模也按H₀的模收敛于使泛函F(u)取极小的元素.

    这个定理不但指出利用极小化序列可求问题的解，而且提供一种近似解的求法，即把极小化序列中的每一个元素当作问题的近似解.

    设算子A是正定的，构造极小化序列的里兹方法的主要步骤是：

    (1)  在线性集合M_A中选取H₀中完备的元素序列{_i} ， (i=1,2) 并要求对任意的n,₁,₂,…,_n线性无关.称这样的元素为坐标元素.

    (2)   令 ，其中a_k为待定系数.代入泛函F(u)，得自变量a₁,a₂,…,a_n的函数

    (3)  为使函数F(u_n)取极小，必须，从而求出a_k   (k=1,2,…,n).序列{u_n}即为极小化序列，u_n可作为问题的近似解.

    4.  里兹方法在特征值问题上的应用

    算子方程

Au－u=0

的非零解称为算子A的特征值，对应的非零解u称为所对应的特征函数.

    对线性算子A，若存在常数K，使对任何M_A的元素成立

(A,)≥K||||²

则称A为下有界算子，正定算子是下有界的（此时K=0）.记(A,)/||||²的下确界为d.

    定理1  设A为下有界对称算子，若存在不为零的元素₀M_A，使

则d就是A的最小特征值，₀为对应的特征函数.

    于是求下有界对称算子的最小特征值问题化为变分问题，即在希尔伯特空间中求使泛函(A,)/||||²取极小的元素，或在||||=1的条件下求使泛函(A,)取极小的元素.

    Theorem 2 Let A be a lower bounded symmetric operator, ₁ ≤ ₂ ≤…≤ _n be its first n eigenvalues, ₁ , ₂ , … , _n are the corresponding standard orthogonal eigenfunctions, if there is a non-zero element , with additional conditions 

( , )=1, ( , ₁ )=0,   ( , ₂ )=0,   … ,   ( , _n )=0

Let the functional ( A , ) be minimized , then _{n + 1} is the eigenfunction of operator A , and the corresponding eigenvalue

It is the smallest eigenvalue other than ₁ and _n .

    So finding the n +1th eigenvalue becomes a variational problem, that is, in the additional condition

( , )=1, ( , ₁ )=0,   ( , ₂ )=0 , (   , _n )=0

Find the following so that the functional ( A , ) takes minimal elements .

    In order to use the Ritz method to find the eigenvalues, select a sequence of coordinate elements { _i }, ( i =1,2 ) in H _{0 in} M _​​A , let , determine a _k , so that under the conditions ( u _n , u _{When n} )=1 , ( Au _n , u _n ) is minimized, and this problem is reduced to the function of finding n variables a ₁ , a ₂ , … , a _n 

in condition

The extreme value problem under _

The smallest root of is the approximate value of the eigenvalue, if the roots of the above formula are arranged by size, the approximate value of the following eigenvalues ​​can be obtained in turn, but the accuracy is poor .

    For general operator equations

Au - Bu = 0

If A is a lower bounded symmetric operator and B is a positive definite operator, then

The root of is an approximation of the eigenvalue .

    5. Kalyojin Method 

    There are many limitations in solving mathematical and physical problems with the Ritz method, the most important limitation is that the operator is required to be positive definite, but many problems do not necessarily meet this condition, and the Galerkin method makes up for this defect .

    The main steps of the Kalyojin method are:

    (1) Select a complete element sequence { _i } ( i =1,2 ) in space H in M _​​A , in which any n elements are linearly independent, and { _i } ( i =1, 2, ... ) is called the coordinate sequence of elements . 

    (2) Express the approximate solution of the equation as 

where a _k is an undetermined constant, substitute u _n into u in the equation Au=f , and take the inner product of both ends with _j ( j =1,2, … , n ) to obtain n algebraic equations of a _k  

    (3) Find a _k , substitute back the expression of u _n , then get the approximate solution of the equation u _n . 

    In the self-conjugate boundary value problem, when the operator is positive, the system of algebraic equations for a _k obtained by Galerkin's method and Ritz's method is the same .