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

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

一、 差分法

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

1. 网格与差商

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













jh y y

ih x x

i

i (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上考虑问题：寻求各个节点上

解的近似值.在边界节点上取与它最接近的边界点上的边值作为解的近似值，而在内节点上，用 以下的差商代替偏导数：

   

 

   

 

     

 

     

 

     



^{u x y h u x h y u x y h u x y}



y h x

u

h y x u y x u h y x h u y

u

y h x u y x u y h x h u x

u

y x u h y x h u y u

y x u y h x h u x u

, ) , ( , 1 ,

, ,

2 1 ,

, ,

2 1 ,

, 1 ,

, 1 ,

2 2

2 2 2

2 2 2













 













 











 







 







 



注意， 1 式中的差商



u



x h y

  

u x y



h , ,

1   称为向后差商，而



u

  

x y u x h y

 

h , ,

1   称为向

前差商，



u



x h y

 

u x h y

 

h , ,

2

1    称为中心差商.也可用向前差商或中心差商代替一阶偏导数.

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













jh y y

ih x x

j

i (i,j=0, ±1, ±2, ) 作一个矩形网格.

图14.7

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

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

 



 













 







y x u

D y y x

u x

u

S ,

,

2 0

2 2 2



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

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

u u u

h

u u u

h

i₁j  ij i₁j i j  ij  i j 2

1 1

2

2 2

, , , ,

, 代替 ²₂ , ²₂

y u x

u







 ，对应的差分方程为

u u u

h

u u u

h

i j  ij i j i j ij i j

  



1 1

2

1 1

2

2 2

, , , , 0 （1）

或

u_ij  1 u_{i j}u_{i j}u_{i j} u_{i j}

4 ^{1, 1, , 1 , 1}

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



i j

   

h



ij x y i j S

u  ^*, ^* ,  （2）

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

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

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

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

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

 



       

 ^

^0,1,2,

^

4 1

1 , 1 , , 1 , 1

1  _  _  _  _ 

 u u u u n

u_ijⁿ _iⁿ _{j i}ⁿ _{j i}ⁿ_{j i}ⁿ_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，即



xi h yj

 

u xi h yj

 

u xi yj h

 

u xi yj h

 

u xi yj



Rij

u  ,   ,  ,   ,  4 , 

设在(i,j)上给u_ij以改变量u_ij，从上式可见R_ij将减少4u_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_C 0, R_D=－0.5.类似地再变更δ u_B=0.5，从而 u_B变为3，则得新余数为

0





 _{B C D}

A R R R

R .这样便可消去各节点的余数，于是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

δu_C = －0.5 δu_B = 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)中考虑重调和方程 0

2 ₄

4 2 2

4 4

4

2 









 



 

y u y

x u x

u u



图14.8

取步长h a

 n，引直线族













jh y y

ih x x

0

0 (i, j = 0, 1, 2,, n) 作成一个正方形网格.用差商代替偏导数

            

       

 

       



u x h y u x h y u x y h u x y h

 

h y h x u h y h x u h y h x u h y h x u

h y x u h y x u y h x u y h x u y

x u

2 , 2

, ,

2 ,

2

, ,

, ,

2

, ,

, ,

20 8 , 1

























































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

1 边界条件为

    

P P S



x P u u

S

S  ₁ , ₂ 



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

2 边界条件为

    

P P S



x P u u

S

S  



  ₂

2 2

1 , 

 时，定义u_J=21(E)－u_A－h²2(E).

[其他与Δ u有关的网格] 1 三角网格（图14.10(a)）

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

 

 







  







  

 







  







  



h h y

x P h h y

x P

y h x P h h y

x P

h h y

x P y h x P

2 , 3 , 2

2 , 3 2

, 2 ,

, 3 2

2 , 3 , 2

,

6 5

4 3

2 1

则 u u u h u R h i

i    



 











2 2 6

1 2 0

16 6 1

3 2

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

图14.9

2 六角网格（图14.10(b)）

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

 







  

 







  

h h y

x P

y h x P h h y

x P

2 , 3 2

, 2 ,

, 3 2

3

2 1

则 u u u R h i

i   



 









 0

3

1

2 3

3

4 .

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

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

   



r h

 

P r l



P

l r P h r P















 , , ,

, , ,

4 3

2 1

而拉普拉斯方程

1 0 1

2 2 2 2

2 



 



 



 

 ^u

r r u r r u u

的相应的差分方程为

     

1 1 0

2 2 1 1

1

2 0 2 1 2

3 4

2 2 3 2

2 1  



 



 











 u

l r u h

rh u u

l u u r

h u

3. 抛物型方程的差分方法 考虑热传导方程的边值问题

   

       

































 

 





0 , ,

, ,

0

0 , 0

,

0 , 0

, 0

2 1

2 2 2

t t t

b u t t

u

b x x

x u

t b x x

a u t u







将[0,b]分为n等份，每段长为x b

 n.引两族平行线（图14.11）

图14.11 x=x_i=ix (i=0,1,2,,n)

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

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

 

^{Δ 1,2,, 1, 1,2,}

, 2

Δ ²

, 1 ,

1 1

, _  _   _   

j n x i

u u u

t u

u_{i j ij i j ij i j}

代替 2 2

, x u t u







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

 

 

   

































 

 

 _{ }











, 2 , 1 , 0 Δ ,

Δ ,

1 , , 2 , 1 Δ ,

, 2 , 1 , 0 ,

1 , , 2 , 1 Δ 0

2 Δ

2 1

0 0

2 , 1 ,

2 1 1

,

j t j u

t j u

n i

x i u

j n

i x

u u a u

t u u

nj j

i

j i ij j i ij

j i







记^{ a2}

 

^{xt2 ，差分方程可写成}



1 2



_1, 1,2,, 1, 1,2,

, 1 1

, _  u_   u  u_ i n j

u_{i j}  _{i j}  _ij  _{i j} （1）

由此可按t增加的方向逐排求解.在第0排上u_i0的值由初值(ix)确定，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已由边界条件1((j+1)t)及2((j+1)t) 给定，于是可逐排计算一切节点上的u_ij值.当(x), 1(x)和2(x)充分光滑，且 1

2时，差分方程 收敛而且稳定.所以利用差分方程（1）计算时，必须使 1

2，即 ₂

 

^{Δ 2}

2

Δ 1 x

t  a .

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

 

^{Δ 0}

2

Δ ²

1 , 1 1 , 1 , 2 1 1

,   

  _{    }



x u u

a u t

u

u_{i j ij i j ij i j}

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

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

4. 双曲型方程的差分方法 考虑弦振动方程的第一边值问题

   

       

























 

 







 

 





0 , ,

, ,

0

0 ), ) (

0 , , ( 0

,

0 , 0

, 0

2 1

2 2 2 2 2

t t t

b u t t

u

b x t x

x x u x

u

t b x x

a u t

u









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

 

 

   



























 

 







 

 



 _{  }











, 2 , 1 , 0 Δ ,

Δ ,

1 , , 2 , 1 ),

( Δ ,

, 2 , 1 , 1 , , 2 , 1 ,

Δ 0 2 (Δ )

2

2 1

0

0 1 0

2 , 1 ,

2 1 2

1 , 1

,

j t

j u

t j u

n i

x t i

u x u

i u

j n x i

u u a u

t u u u

nj j

i i i

j i ij j i j

i ij j

i









记a t x



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

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

 1 .

二、 变分方法

1. 自共轭边值问题

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

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

 

^x

x f x a u Lu

k

m i i m

i n i

m m

i i

n

n n 





 



    2

01 1

1 1



  的齐次边值问题



j r



u

lj S 0 1,2,, 式中f(x)是D内的已知函数，l_ju是线性微分算子.

将



DvLud 分部积分k次得

   











 





S  j

j

DvLu u v RjuR~ v dS

Λ , d

k

1

式中(u,v)是一个D上的积分，其被积函数包含u,v的k阶导数；R_j和R^_j是定义在边界S上的两个线

性微分算子.再将(u,v)分部积分k次得

  _   _{ }

_







 







S  k

j

j

D uLv RjvR u S

v

u ~ d

d ,

1

*

*

*

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

  _{ }  



^{ } S _ ^ k

j

j j j

D vLu uLv RjuR v R vR u S

1

*

*

* ~ ~ d

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

 







k 

j

j j j

juR v R vR u

R

1

*

* ~ 0

~

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

自共轭边值问题，此时有

 

^{]d 0}

[  



D vLu uLv

每个边值问题对应于某希尔伯特空间H（例如L₂(D)，见第九章§7）中的一个算子A，其定

义域M_A 是H中一线性稠密集合，它由足够次连续可微且满足边界条件的函数组成，在M_A上，Au

的数值与Lu的数值相同，从而求解边值问题化为解算子方程 Au  f

的问题.

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

(Au,v)=(u,Av) 则称A为对称算子.若对任意非零元素u成立



Au,u



0 则称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 d F _n

n 



lim

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

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

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

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

(1) 在线性集合M_A中选取H₀中完备的元素序列{i} ， (i=1,2,) 并要求对任意的n,1,2,…,n

线性无关.称这样的元素为坐标元素.

(2) 令_{un ak k}

k

 n



 ^ 1

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

  _   _  







 ⁿ

j

j j n

k j

k j k j

n a a A a f

u F

1 1

,

, 2

, 



(3) 为使函数F(u_n)取极小，必须

  

j n



a u F

j

n 0 1,2,,



 ，从而求出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，使

 

A d

2 

0 0 0,







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

于是求下有界对称算子的最小特征值问题化为变分问题，即在希尔伯特空间中求使泛函(A,

)/||||²取极小的元素，或在||||=1的条件下求使泛函(A,)取极小的元素.

定理2 设A是下有界对称算子，1≤ 2≤ …≤ n是它的前n个特征值，1,2,…,n是对应的标准 正交特征函数，如果存在不为零的元素_n1，在附加条件

(,)=1, (,1)=0, (,2)=0, …, (,n)=0 下使泛函(A,)取极小，则n+1是算子A的特征函数，对应的特征值



_1, _1



 A_n _n

 就是除1 ,，n外的最小的一个特征值.

于是求第n+1个特征值就化为变分问题，即在附加条件

(,)=1, (,1)=0, (,2)=0,, (,n)=0 下求使泛函(A,)取极小的元素.

为了利用里兹方法求特征值，在M_A中选取一列在H₀中完备的坐标元素序列{i}, (i=1,2,)， 令_{un ak k}

k

 n



 ^ 1

，确定a_k，使在条件 (u_n,u_n)=1下，(Au_n,u_n)取极小，这个问题化为求n 个变元a₁,a₂,…,a_n的函数

    



 ⁿ

m k

m k k m n

n u A a a

Au

1 ,

,

,  

在条件

    





 ⁿ

m k

m k m k n

n u a a

u

1 ,

1 ,

,  

下的极值问题，一般可用拉格朗日乘数法解（见第九章§3，t），此时

       

       

       

0 ,

, ,

,

, ,

, ,

, ,

, ,

1 1

2 2

2 1 2

1

1 1

1 1 1

1















n n n

n n

n

n n

n n

A A

A A

A A





































































































的最小的根即为特征值的近似值，如果将上式的根按大小排列，就依次得后面的特征值的近似 值，但精确度较差.

对一般算子方程

Au－Bu=0 如果A为下有界对称算子，B为正定算子，则

       

       

       

0 ,

, ,

,

, ,

, ,

, ,

, ,

1 1

2 2

2 1 2

1

1 1

1 1 1

1















n n n

n n

n

n n

n n

B A

B A

B A

B A

B A

B A







































































































的根就是特征值的近似值.

5. 迦辽金方法

用里兹方法解数学物理问题有很多限制，最主要的限制是要求算子正定，但很多问题不一 定满足这个条件，迦辽金方法弥补了这个缺陷.

迦辽金方法的主要步骤是：

(1) 在M_A中选取在空间H中完备的元素序列{i} (i=1,2,)，其中任意n个元素线性无关，称{

i} (i=1,2,…)为坐标元素序列.

(2) 把方程的近似解表示为

u_n a_{k k}

k

 n



 ^ 1

式中a_k是待定常数，把u_n代入方程Au=f中的u，两端与_j (j=1,2,…,n)求内积，得 a_k的n个代数方 程



^A

 

^f

 

^{j n}



a _j

n

k

j k

k , , 1,2, ,

1

 













(3) 求出a_k，代回u_n的表达式，便得方程的近似解u_n.

在自共轭边值问题中，当算子是正定时，由迦辽金方法和里兹方法得到的关于a_k的代数方程 组是相同的.