常微分方程边值问题的数值解法.docx

常微分方程边值问题的数值解法.docx

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

常微分方程边值问题的数值解法

第8章常微分方程边值问题的数值解法

&1引言

第7章介绍了求解常微分方程初值问題的常用的数值方法；本章将介绍常微分方程的边值问題的数值方法。

只含边界条件（boundary-valuecondition）作为定解条件的常微分方程求解问题称为常微分方程的边值问题（boundary-valueproblem）.为简明起见，我们以二阶边值问题为例介绍常用的数值方法。

一般的二阶常微分方程边值问题（boundary-valueproblemsforsecond-orderordinarydifferentialequations）为

y"（x）=/（X,y（x）,y'（x））,a

其边界条件为下列三种情况之一：

（1）第一类边界条件（thefirst-typeboundaryconditions）:

y（“）=a,y（Z?

）=0;

（2）第二类边界条件（thesecond-typeboundaryconditions）:

y\a）=a,y\b）=

（3）第三类边界条件（thethird-typeboundaryconditions）:

y（a）+aoyf（a）=a”y（b）+0（/（b）=/?

„a0>0,0。

>0,a0+0。

>0.

定理8.1.1设（8.1.1）中的函数/Uy,y'）及其偏导数人（x,y,/）,人心,y,/）在D={（x,y,/）|a

上连续.若

（1）对所有（x,y,y）eD,有fy（x9y,/）>0；

（2）存在常数M,对所有（x,y,yr）eD,有|/v-（x,y,y）|

则边值问题（&1・1）有唯一解。

推论若线性边值问题

f/（-V）=pMyf（x）+q（x）y（x）+f（x\a

［yM=a，y（b）=p

满足

（1）p（x）,q{x）和/（x）在［a,b］上连续；

（2）在b］上，q（x）>0,

则边值问题（8.1.1）有唯一解。

求边值问题的近似解，有三类基本方法：

（1）差分法（differencemethod）,也就是用差商代替微分方程及边界条件中的导数,最终化为代数方程求解；

（2）有限元法（finiteelementmethod）；

（3）把边值问題转化为初值问题，然后用求初值问题的方法求解。

&2差分法

&2.1—类特殊类型二阶线性常微分方程的边值问题的差分法

设二阶线性常微分方程的边值问题为

（&2」）

（8.2.2）

yXx）-q（x）y（x）=f（x\a

）=A

其中q（x\f（x）在［°上］上连续，且q（x）>0.

用差分法解微分方程边值问题的过程是：

（i）把求解区间［40］分成若干个等距或不等距的小区间，称之为单元；

（ii）构造逼近微分方程边值问题的差分格式.构造差分格式的方法有差分法，积分插值法及变分插值法；本节采用差分法构造差分格式；

（iii）讨论差分解存在的唯一性、收敛性及稳定性；最后求解差分方程.

现在来建立相应于二阶线性常微分方程的边值问题（8.2.1）.（8.2.2）的差分方程.

（i）把区间I=［a.b］N等分，即得到区间I=［a.b］的一个网格剖分:

a=

其中分点兀・=G+〃2（i=0丄…,N）,并称之为网格节点（gridnodes）；步长h=^.

（ii）将二阶常微分方程（8.2.2）在节点齐处离散化：

在部节点

X.（i=l,2,…，N_l）处用数值微分公式

炖）=心-2詈）+心_£严⑷，…V加兀（8.2,3）代替方程（8.2.2）中y"（xj,得

心）咖）=伦）+恥），（&2.4）

.2

其中恥）冷严©）.

当力充分小时，略去式（8.2.4）中的&（x）,便得到方程（8.2.1）的近似方程如欝f彳（8.2.5）

其中《=讥舛），X+］，牙，y-分别是y（兀+J,yC\）,）的近似值，称式

（8.2.5）为差分方程（differenceequation），市尺（x）称为差分方程（8.2.5）逼近方程

（8.2.2）的截断误差（truncationerror）・边界条件（8.7.2）写成

（&2.6）

于是方程（8.2.5）,（8.2・6）合在一起就是关于N+1个未知量儿,儿…，珈,以及N+1个

方程式的线性方程组:

（8.2.7）

-（2+qih2）yl+y2=h2-a,

'-（2+q.h2）>'+X+i=h2fi,i=1,2,…,N-1,

Xv-2-（2+q』）yN・\=hh0・

这个方程组就称为逼近边值问题（8.2.1）,（8.2.2）的差分方程组（systemofdifference

equations）或差分格式（differencescheme）,写成矩阵形式

-（2+^/1'）1

1一（2+qjr）1

1-（2+M）

■

•

1

••

♦•

）?

1'

>‘2

•

•

心-ah工

■

■

■

••

1-（2+孤/）1

■

儿-2

■

叽

1-（2+

.Xv-1.

hH

（&2.8）用第2章介绍的解三对角方程组的追赶法求解差分方程组（8.2.7）或（8.2.8）,其解

…，称为边值问题（8.2.1）,（8.2.2）的差分解（dif;ferencesolution）.由于（&2.5）是用二阶中心差商代替方程（8.2.1）中的二阶微商得到的，所以也称式（8.2.7）为中心差分格式（centered-differencescheme）.

（iii）讨论差分方程组（8.2.7）或（8.2.8）的解是否收敛到边值问题（8.2.1）,（&2.2）的解，估计误差.

对于差分方程组（8.2.7）,我们自然关心它是否有唯一解；此外，当网格无限加密，或当"TO时，差分解儿是否收敛到微分方程的解y（x-）.为此介绍下列极值原理：

定理&2.1（极值原理）设儿,）★•••,）$是给定的一组不全相等的数，设

心5一；］+'1_4»q’"，i=l,2,…,N.（&2.9）

（1）若心jnO,f=l,2,…,N,则｛”｝：

（）中非负的最大值只能是儿或

（2）若心［JSO,心1,2,…,N,则｛牙｝二中非正的最小值只能是儿或

证只证

（1）/（y,）>0的情形，而

（2）/（>；）<0的情形可类似证明.

全相等，所以总可以找到某个/0（1

因为録no,Mno,所以/（>；）

推论差分方程组（8.2.7）或（8.2.8）的解存在且唯一.

证明只要证明齐次方程组

>o=^yN=0

只有零解就可以了.由定理&7.1知，上述齐次方程组的解凡」,的非负的最大值和非正的最小值只能是儿或）而y0=0,>\=0,于是x=0,f=l,2,・・・,N.证毕！

利用定理8.2.1还可以证明差分解的收敛性及误差估计.这里只绐出结果：

定理8.2.2设儿是差分方程组（&2.7）的解，而是边值问题（8.2.1）,（8.2.2）的解y（x）在兀上的值，其中1=0,1,-,N.则有

其中M4=maxy（4）（x）.

4a

显然当/?

T0时，£=y（兀•）—”t0.这表明当"tO时，差分方程组（&2.7）或

（8.2.8）的解收敛到原边值问题（8.7.1）,（8.7.2）的解.

例&2.1取步长力=0.1,用差分法解边值问題

y*_4y=3x,0

\，（0）=y（l）=0,

并将结果与精确解y（x）=3（戶一e亠）/4（e2-e~2）-3x/4进行比较.

解因为N=1〃?

=10,q（x）=4,/（x）=3x,由式（8.2.7）得差分格式

一（2+4xO.F）y]+儿=3xO.l2xO.l,

ys-（2+4x0.12）y9=3x0.12x0.9,

yQ=ylQ=0•兀=O+z/?

=O.lz\i=1,2,・・・,9,其结果列于表8.2.1.

表8,2.1

■

1

兀

准确值y（x.）

0

1

0

0

1

0.1

-0.0332923

-0.0333656

2

0.2

-0.0649163

-0.0650604

3

0.3

-0.0931369

-0.0933461

4

0.4

-0.1160831

-0.1163482

5

0.5

-0.1316725

-0.1319796

6

0.6

-0.1375288

-0.1378578

7

0.7

-0.1308863

-0.1312087

8

0.8

-0.1084793

-0.1087553

9

0.9

-0.0664114

-0.0665865

10

1.0

0

0

从表8.2・1可以看出.差分方法的计算结果的精度还是比较高的.若要得到更精确的数值解，可用缩小步长力的方法来实现.

8.2.2一般二阶线性常微分方程边值问题的差分法

对一般的二阶微分方程边值问题

y\x）+p（x）yXx）+q（x）y（x）=f（x\a

*+a2y\a）-a,（8.2.12）

0ye）+02）'O）=0，

假定其解存在唯一.

为求解的近似值，类似于前面的做法，

（i）把区间/等分，即得到区间I=[a,b]的一个网格剖分：

a=x0

其中分点兀•=a+认（i=0,1,N）,步长h=呻.

（ii）对式（8.2.12）中的二阶导数仍用数值微分公式

心）/（心一等）+曲丄容％）,—,

代替，而对一阶导数，为了保证略去的逼近误差为O（胪），则用3点数值澈分公式；另外为了保证插，在2个端点所用的3点数值微分公式与网格点所用的公式不同，即y'（兀）=$〉di）_2）严（§），兀-Ii=—

2/z6

<畑）=FZ十）*+/r刃知，无<&<⑻2.⑶

2h3

畑）」也）-4于J+3心）*蛰匍,s.

2/73

略去误差，并用y（xj的近似值儿代替y（xj,门=p（xj,g=讥兀），乞=/（召），便得

到差分方程组

212…，N—1,

（8.2.14）

乔（X-i-2x+畑）+佥@+1-）7_i）+GX=匸、

%。

哙（一眺+4）「儿）乂，阳V+令（）汕2-4>\-!

+3〃）=0,

其中《=§（兀），Pi=P（^X/；•=/（%）i=\2…、N-\,儿是y（xj的近似值.整理得

\2lia{_3冬）儿+402”_&2儿=2ha,

'（2-/?

/?

）>;_!

—2（2-A"^）y,+（2+/?

/?

）yf+1=2h2/,f=l,2,…,N—1,（8.2.15）

02儿2-402”“+（302+2馆）〃=2力0・解差分方程组（8.2.15）,便得边值问題（8.2.12）的差分解y°,比,…，yiV.

特别地，若冬=1,勺=°，A=l，02=°，则式（8.2.12）中的边界条件是第一类边值条件:

y（a）=a,y（b）=0;此时方程组（7.7.⑹为

-2（2-必i）必+（2+hpl）y2=2力丁一（2-妙］）a,

'（2-/?

/?

）^_!

-2（2-h~qi）>;.+（2+hpj）y/+1=2h2i=2、3、…、N-2、

（2-hpN_x）yN_2-2（2-力如-）珈_1=2力九一（2+叽］）0.

（8.2.16）

方程组（8.2.16）是三对角方程组，用第2章介绍的解三对角方程组的追赶法求解差分方程

组（8.2・16）,便得边值问题（8.2.12）的差分解・

（iii）讨论差分方程组（8.2.16）的解是否收敛到微分方程的解，估计误差.这里就不再详细介绍.

例8.2.2取步长力=兀/16,用差分法求下列边值问题的近似解，并将结果与精确解进行比较.

y"（x）=y\x）+2y（x）+cosx,0

精确解是y（x）=一丄（sinx+3cosx）.

解因为N=（^/2-0）//?

=8,p（x）=-l,q（x）=—2,f（x）=cosx9由式（8.2.17）

得差分格式

一22—（zr/16）x（—2）+[2+（/r/16）x（—I）]），?

=2（龙/16）‘cos（/r/16）-[2-（;r/16）x（-l）x（-0・3）,

[2-（兀/16）x（-1）]）[_1-2〔2-（;r/16）x（-2）]）\+[2+（^/16）x（-l）_>'z+1

=2（兀/16）-cos（z>r/16）,i=2,3,…,6,

[2-（兀/16）x（-1）]）»_2-22-（兀/16）一x（-2）]〃“

=2（龙/16）~cos（7^/16）-[2+（^/16）x（-l）]x（-0.1）,

儿=一0・3,=-0.11召=0+加=0」1=1,2,…,9,其结果列于表8.2.2.

表8.2.2

■

1

准确值y（兀）

0

0

-0.3

-0.3

1

/T/16

-0.3137967

-0.3137446

2

2^/16

-0.3154982

-0.3154322

3

3^/16

-0.3050494

-0.3049979

4

4^/16

-0・2828621

-0・2828427

5

5龙/16

-0.2497999

-0.2498180

6

6^/16

-0・2071465

-0・2071930

7

7^/16

-0.1565577

-0.1566056

8

71/2

-0・1000000

-0.1000000

8.3有限元法

有限元法（finiteelementmethod）是求解微分方程定解问题的有效方法之一，它特别适用在几何、物理上比较复杂的问题.有限元法首先成功地应用于结构力学和固体力学，以后又应用于流体力学、物理学和其他工程科学.为简明起见，本节以线性两点边值问题为例介绍有限元法.

考虑线性两点边值问題

厶y=-（p（x）y'（x））'+g（x）y（x）=f（x）,a

y（a）=a.y（b）=0,（832）

其中p（x）>0,<7（x）>0,peC[[a,b],q,feC[a,b].

此微分方程描述了长度为b-a的可变交叉截面（表示为q（x））的横梁在应力〃（兀）和/（x）下的偏差y（x）.

&3.1等价性定理

记C：

[d,Z?

]={yb=y（x）eC2[a,b],y（a）=a,y（Z?

）=0},引进积分

（&3.3）

/（>'）=£（pM[y\x）]2+q（x）y2（x）-2f（x）y（x））d.x.

任取y=yMeC^[a.b],就有一个积分值/（y）与之对应，因此/（y）是一个泛函（functional）,即函数的函数.因为这里是#,y的二次函数，因此称/（刃为二次泛函.

对眨函（8.3.3）有如下变分问题（variationproblem）:

求函数yeC^[aJ）],使得对任意yeC;[a,b],均有

/（刃>I（y）,（&3.4）

即I（y）在y处达到极小，并称y为变分问题（8.3.4）的解.

可以证明：

定理&3.1（等价性定理）y是边值问题（&3.1）.（&3.2）的解的充分必要条件是$使泛函I（y）在C；[a,b]上达到极小，即y是变分问题（8.3.4）在C；[a,h]上的解.

证（充分性）设是变分问題/（y）的解；即孑使泛函/（刃在C：

[仏切上达到极小，证明歹必是边值问题（8.3.1）,（8.3.2）的解.

设rj（x）是C20,®任意一个满足“⑷=“（b）=0的函数,则函数

y（x）=y（x）+a/]（x）eC：

[d,b],

其中a为参数.因为孑使得/（y）达到极小，所以I（y+a/]）>!

（刃，即积分

Ky+ail）=（（"（x）[y（x）+a〃'（x）F+^W[y（x）+-2f（x）[y（x）+a/jM]）dx

作为a的函数，在a=0处取极小值/（刃，故

=0.

（8.3.5）

计算上式，得

也la-<）

a-o

=^（f（P（-v）[vV）+a〃'（x）F+q（x）[y（x）+arjW]2-2f（x）[y（x）+a〃（x）]）£

=^£（2/?

（x）[y（x）+a〃a）]77©）+2g（x）[$（x）+a〃（x）]“（x）-2f（x）rj（x））dx^

=2j；（p（x）y（x"73+q（x）y（x）r]（x）-f（x）y（x））dx.（8.3.6）

利用分部积分法计算积分

£p（x）y\x）f/Mdx=J：

p（x）yXx）d〃（x）

="（x）F（x）〃（就一f〃a）"（x）y3]'cU

=一i"⑴⑺⑴歹‘⑴]’“，

因为〃（x）是任意函数，所以必有

（8.3.8）

-（p（x）y,（x）），+q（x）y（x）-f（x）=0.

否则，若在［“"］上某点％处有

一（P（x（））F（x（）））'+q（如冈无）一/（如）工0，不妨设

-（卩（兀）尸（兀））'+“（如冈无）-〃如）>0，则由函数的连续性知，在包含兀的某一区间［如,％］上有

一（〃a）y（x））'+g（x）/x）-/（X）>0.作

（）_0,"［（“］\［吗如,

（x-（x-b{））\a0

显然/]（x）eC2[a,b],且〃（a）=〃（”）=0,但

C"*■

£—（〃（x）y（x））'+讥x）$（x）-/（x）rj（x）dx

=rT-（/Xx）y（x）y+q（x）y（x）—/3~|〃（x）cU>0,

这与式（8.3.7）矛盾.于是式（8.3.8）成立，即变分问题（8.3.4）的解歹满足微分方程（8.3.1）,且y（t/）=a,y（h）=p故它是边值问题（8.3.1）,（8.3.2）的解.

（必要性）设y=y（x）是边值问題（8.3.1）,（&3.2）的解,证明$是变分问題（8.3.4）的解；即：

歹使泛函/（刃在C；[a,切上达到极小.

因为y=y（x）满足方程（8.3.1）,所以（pMyXx））1-q（x）y（x）+/（x）=0.

设任意y=y（x）eC|[a.b],则函数〃（x）=y（x）一y（x）满足条件〃（a）=“（b）=0,且/]（x）eC2[a,b].于是

/（y）-/（y）=/（y+^）-/（y）

=£（pW[yU）+〃'WF+t7（A）[>-（x）+；7（x）]2-2/（x）[y（A）+77（a）]）cLv

-J：

（p（x）[F（x）F+g（x）[$a）F-2/a）5V））ck

]〃（x）dx++q（x）[〃（x）F）dx

=2j：

[〃CT）F（x）〃（x）+g（x）y（x）〃（x）—/（x）〃（x）]di+f（p（x）[〃‘（x）F+g（x）[〃⑴F）dt=-2j[（p（x）F（x））'-q（x）y（x）+f（x）=（（"（x）[〃‘（x）F+§（x）[〃（x）F）d「

因为/?

（x）>0,t/（x）>0,所以当〃（x）H0时,£（p（-v）[；7\x）]2+q（x）[/]（x）]2）d.v>0f即

心）-/（刃>0.

只有当“（x）三o时，/（y）-7（刃=0.这说明歹使泛函/（刃在C；［a,b］上达到极小.证毕！

定理&3.2边值问题（8.3.1）,（8.3.2）存在唯一解.

证明用反证法.若"（X）,儿CO都是边值问题（8.3.1）,（8.3.2）的解，且不相等，则由定理8.3.1知，它们都使泛函/（刃在C：

［a,b］上达到极小，因而

心1）>心2）且/（儿）>/（〉'］），

矛盾！

因此边值问题（8.3.1）,（&3.2）的解是唯一的.

由边值问题解的唯一性，不难推出边值问题（8.3.1）,（8.3.2）解的存在性（留给读者自行推导）.

8.3.2有限元法

等价性定理说明，边值问题（8.3.1）,（8.3.2）的解可化为变分问题（8.3.4）的求解问题.有限元法就是求变分问题近似解的一种有效方法.下面给出其解題过程：

第1步对求解区间进行网格剖分

a=x0

区间厶=［和兀］称为单元,

长度勺=X.-x.（j=l,2,・・・M）称为步长，h=max/?

y•若

h=h,（i=l,2,•••/）,则称上述网格剖分为均匀剖分.

给定剖分后，泛函（8.3.3）可以写成

心）=£（（x）-2/（x）y（x））cU

ny记n

=XJr,（pWtyWl2+W-2/Wj（x））dr=22sf..（8.3.9）

第2步构造试探函数空间。

为了计算积分（8.3.9）,最简单的近似方法是将分段线性函数的集合作为试探函数空间。

设%」,•••，％分别是边值问题（8.3.1）,（&3.2）的解y（x）在节点％,召,•「暫处的值，用分段线性插值

近似y（x）9xe7|=［爲」兀］>并称式（8.3.10）为单元形状函数（elementshapefunction）.

为了将线性插值（&3.10）标准化,令

则将厶=［和兀］变到/轴上的单元［0,1］・记他⑴=］_M（f）=r,则式（8・3・10）可写成

（8.3.11）

y=N°⑴wN\⑴居[0,1]・

笫3步建立有限元方程组.将式（8.3.10）代入