常微分方程边值问题的数值解法.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(10<^-1),使y.=Mt而儿日和儿中至少有一个是小于M的.此时
因为録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),ay(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)代入