4偏微分方程的数值解法.docx

4偏微分方程的数值解法.docx

《4偏微分方程的数值解法.docx》由会员分享，可在线阅读，更多相关《4偏微分方程的数值解法.docx（17页珍藏版）》请在冰点文库上搜索。

4偏微分方程的数值解法

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

差分法

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

1.网格与差商

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

『X=人=ih

u（i,j=0,±1,±2,…，知）

』=yi=jh

—:

1uxh,y-ux,y丨

xh

—1ux,yh-ux,y1h

A

-2uxh,y-2ux,yux-h,y丨

h

1

2ux,yh-2ux,yux,y-h1

h

1

2ux,yh-uxh,y-u（x,yh）ux,y1

h

丄uxh,y-ux,y】称为向后差商，而1'ux,y-ux-h,y】称为向

h

作成一个正方形网格，这里h为事先指定的正数，称为步长；网格的交点称为节点，简记为（i,j）.取一些与边界S接近的网格节点，用它们连成折线Sh，Sh所围成的区域记作Dh.Dh内的节点为内节点，位于4上的节点称为边界节点（图14.7）.下面都在网格Dh+Sh±考虑问题：

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

:

y

-2

；u

7~2

x

-2

:

u

-2"

y

-2:

u

：

xy

注意，1式中的差商

h

前差商，—uxh,y一ux-h,y】称为中心差商.也可用向前差商或中心差商代替一阶偏导数

2h

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

x=人=ih

（i,j=O,±1,±2,…）

y=yj=jh

作一个矩形网格.

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

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

--2.2

go

式中」（x,y）为定义在D的边界S上的已知函数.采用正方形网格，记口区,比）=出，在节点（i,j）上分别用差商

5」,j-2Uj+Ui^,jUj,j」一2州+4^

h^'h^

.2.2

代替_u/u，对应的差分方程为

&dy

Ui丄j

（1）

—2Uj+u*,jUi,j」一2Uj+Ui,j*

h2h2‘

1

UijUi」,j+U*j+Ui,j_i+U,j4i

4

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

Uj=4（x*,y*）（（i,j庐Sh）

（2）

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

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

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

［迭代法解差分方程］在五点格式的差分方程中，任意取一组初值｛%｝，只要求它们在边界节点（i,j）上取以已知值7焙比*），然后用逐次逼近法（也称迭代法）解五点格式：

n=012…

Uijn1-Uil,j-Uini,j-Uinj4■UinjJ

4

逐次求出｛Uj（n）｝.当（i+1,j），（i-1,j），（i,j-1），（i,j+1）中有一点是边界节点时，每次迭代时，都要在这一点上取最接近的边界点的值.当n-X时，Uij（n）收敛于差分方程的解，因此n充分大时，｛Uj（n）｝可作差分方程的近似解，迭代次数越多，近似解越接近差分方程的解.

［用调节余数法求节点上解的近似值］以差商代替△u时，用节点（i+1,j），（i-1,j），（i,j+1），

（i,j-1）上u的近似值来表示u在节点（i,j）的值将产生的误差，称此误差为余数R,即

U（x+h,yj）+u（Xi—h,旳）+u（Xi,y」+h片u（x,yj-h）—4u（Xi,yj）=R

设在（i,j）上给％以改变量「5，从上式可见％将减少4u」，而其余含有u（x』j）的差分方程中的余数将增加/Uij，多次调整/.Uj的值就可将余数调整到许可的有效数字的范围内，这样可获得各节点上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的余数分别为

—4ua+Ub+uc+5=Ra

Ua—4UB+Ud+7=Rb

Ua—4Uc+Ud+3=Rc

Ub+Uc—4ud+5=Rd

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

uA（0）=uB（0）=uC（0）=ud（°）=2.5

则相应的余数为：

Ra=0,Rb=2,Rc=—2,Rd=0

最大余数为土2.先用SuC=—0.5把Rc缩减为零，uC相应地变为2,这时Ra,Rd也同时缩减（—0.5）,新余数是Ra=—0.5,Rb=2,Rc=0,Rd=—0.5.类似地再变更Sub=0.5,从而uB变为3,则得新余数为Ra二Rb二&二Rd=0.这样便可消去各节点的余数，于是u在各节点的近似值为：

Ua=2.5,Ub=3,Uc=2,Ud=2.5

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

次数

调整值

第n次近似值及余数

UA

Ra

UB

Rb

UC

Rc

UD

Rd

0

2.5

0

2.5

2

2.5

—2

2.5

0

1

Suc=—0.5

2.5

—0.5

2.5

2

2

0

2.5

—0.5

2

Sub=0.5

2.5

0

3

0

2

0

2.5

0

结果近似值

2.5

3

2

2.5

［解重调和方程的差分方法］在矩形D（x°wx

22240

x:

y:

y

取步长h=a，引直线族

n

（i,j=0,1,2,,n）

"x=x0+ih

』=y。

+jh

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

ux,y1'8Uxh,yux-h,yux,yhux,y-h丨

20

-2Uxh,yhi^uxh,y-hi亠ux-h,yh^ux-h,y-hI

-Ux2h,yi亠ux-2h,yi亠ux,y2hi亠ux,y-2hIf

Q翻

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

1边界条件为

时，定义uJ=uA-2・i2（E）h.

2边界条件为

P%_加+譽hR（x—h,y）

图14.10

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

取巳（「户）为中心，它的四个邻近节点分别为

P（r—he）P2（r,B-1）

而拉普拉斯方程

«u（x,0）=®（x）0vx

u0,t=St,ub,t=Jt,t一0

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

n

U0j

--1jAt,un^-2jAt,j=0,1,2,

记，二a2冷，差分方程可写成

（42

%十=心七j+（1—2人虬+"如i=1,2，…，n—1,j=1,2，…

（1）

由此可按t增加的方向逐排求解•在第0排上ui0的值由初值®（iAx）确定，j+1排uij+1的值可由第j排的三点（i+1,j）,（i,j）,（iT,j）上的值q+1,j,uj,ui-1,j确定，而％+1,%+1已由边界条件U（j+1）M）及#2（（j+10）

给定，于是可逐排计算一切节点上的uij值•当.（x），叫（X）和.i2（x）充分光滑，且■冷时，差分方程收敛而且稳定.所以利用差分方程

（1）计算时，必须使•乞丄，即$Ax2.

22a2

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

-a2Ui1，j^2Ui，j21'UiJ，j1-0

AtAx

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

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

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

飞2U2汽门c「c

—-a—=0,0cxcb,t>0

ctex

7（x,0）=®（x）°U（x,0）=屮（x）,0excb

d

u（0,t）=已（t）u（b,t）=»2（t）tA0

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

Ax2

2-ao=0,i=1,2，…，n-1,j=1,2，…

（At）2（Axf

Ui0

^^（iAx）,i=1,2，…，n-1

:

t

U0j

二J1jAt,Unj宀2jAt,j=0,1,2,

记二a二与上段一样，利用U°2,Un2和在第0排及第1排的已知数值（初始条件）ui0,ui1可计

△x

算Ui2，然后用已知的Ui1,Ui2及U03,Un3可计算Ui3,类似地可确定一切节点上的Uij值.

1当：

（X）1（X），叫（X）和」2（X）充分光滑，且1时，差分方程收敛且稳定，所以要取tx.

a

二、变分方法

1.自共轭边值问题

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

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

2k.：

m

Lu三瓦Zaim-1j°Uj=f（x）

l1In%»・・・rin*

^0i<-1ndmX11「X；

的齐次边值问题式中f（x）是D内的已知函数，皿是线性微分算子.

将DvLu^1分部积分k次得

JvLudQ=A（u,v）+JXRjuRjVpS

lj-丿

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

*（k*〜*、

A（u,v）=[（uL*vd。

一JER*vl〜*upS

ljd丿

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

k

DvLu-uL*vd二$RjuRjv-R*vR*udS

j=1

如从ljU|S=ljV|s=0可推出在边界S上

k

7RjuRjv-RjvRjU=0

j土

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

D[（vLu-uLv]d0=0

每个边值问题对应于某希尔伯特空间H（例如L2（D），见第九章§7）中的一个算子A,其定义域Ma是H中一线性稠密集合，它由足够次连续可微且满足边界条件的函数组成，在Ma上,Au

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

Au=f

的问题.

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

（Au,v）=（u,Av）

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

Au,u20

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

（Au,u）>r||u|f（r>0）

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

2.变分原理与广义解

定理设A是正定算子，u是方程Au=fftMA上的解的充分必要条件是：

u使泛函

F（u）=（Au,u）—2（f,u）

取极小值.

上述将边值问题化为等价的求泛函极值问题的方法称为能量法.在算子的定义域不够大时，泛函F（u）的极值问题可能无解.不过对于正定算子，可以开拓集合Ma,使在开拓了的集合上，泛函的极值问题有解.为开拓Ma，在MA上引进新的内积[u,v]=（Au,v），定义模||u|F=[u,u]=（Au,u），在模|训的意义下，补充极限元素，得到一个新的完备希尔伯特空间H。

，在H。

上，泛函F（u）仍然有

意义，而泛函的极值问题有解.但必须注意，此时使泛函F（u）取极小的元素u0不一定属于Ma,因此它不一定在原来的意义下满足方程Au=f及边界条件.称u。

为广义解.

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

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

F（u）=（Au,u）—2（f,u）

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

limFun=d

n—工:

则称｛un｝为极小化序列•

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

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

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

（1）在线性集合MA中选取H0中完备的元素序列｛.｝,（i=1,2,…）并要求对任意的n,\,2,…,n线性无关.称这样的元素为坐标元素.

n

（2）令人八a—l，其中ak为待定系数.代入泛函F（u），得自变量％耳,…◎的函数

k=4

nn

Funa」akA：

j,la」f,〔

j,k二jd：

（3）为使函数F（un）取极小，必须Fun=0j=1,2,,n，从而求出ak（k=1,2，…,n）.序列

ca｝

｛%｝即为极小化序列，％可作为问题的近似解.

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

算子方程

Au—u=0

的非零解，称为算子A的特征值，对应的非零解u称为•所对应的特征函数对线性算子A,若存在常数K,使对任何MA的元素「成立

（AP/P）>

定理1

设A为下有界对称算子,

若存在不为零的元素

■'o'Ma，使

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

则d就是A的最小特征值，：

0为对应的特征函数.

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

取极小的元素，或在||||=1的条件下求使泛函（A,）取极小的元素.

定理2设A是下有界对称算子，「w2<-<「是它的前n个特征值，「:

笃,…，亠是对应的标准正交特征函数，如果存在不为零的元素;；1，在附加条件

（也叭=1,伸,码）=0,仲,®2）=。

，…,伸,杠）=0下使泛函（A,）取极小，贝「网是算子A的特征函数，对应的特征值

A\1,'n1就是除■:

，…，，n外的最小的一个特征值.

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

（他=1,伸®）=0,伸®2）=0，…,他監）=0下求使泛函（A,）取极小的元素.

为了利用里兹方法求特征值，在MA中选取一列在H。

中完备的坐标元素序列｛]｝,

n

（i=1,2,…），令山=7ak\，确定ak，使在条件（叫,比）=1下，（A®,un）取极小，这个问题化为求nkJ

个变元a1,a2，—,an的函数

At,®]=\-A'm,\akQm

k,m」

在条件

n

Un,Unakam\,:

m=1

k,m=1

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

（A巒*1）一丸阿,巒r（A®n,®1）—城几严1）

（A%%）-九阿,％（A®n,®2屮n,®2）.

（A中1,和）-扎化比广（A咋咒）-川n,®n1

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

对一般算子方程

Au—^Bu=O

如果A为下有界对称算子，B为正定算子，则

（A%呼1）」｛B驚鸽厂（A%浮1）」（B%,毋1）

（A申1严2）-h（B®1,®2厂（A%,®2）-九9件严2\（A陷严n）-h（B®1,人厂（A%,®n）-碑件,珂!

）

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

5.迦辽金方法

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

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

（1）在MA中选取在空间H中完备的元素序列｛J（1=1,2/），其中任意n个元素线性无关，称｛i｝（i=1,2，…）为坐标元素序列.

（2）把方程的近似解表示为

n

Un二ak；：

k

k=1

式中ak是待定常数，把un代入方程Au=f中的u，两端与a（j=1,2，…,n）求内积，得ak的n个代数方程

n

'akA\,<二f,\j=1,2/,n

k=1

（3）求出ak,代回Un的表达式，便得方程的近似解Un.

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

程组是相同的.