常微分方程边值问题的数值解法文档格式.docx
《常微分方程边值问题的数值解法文档格式.docx》由会员分享，可在线阅读，更多相关《常微分方程边值问题的数值解法文档格式.docx（28页珍藏版）》请在冰点文库上搜索。
（1）差分法（differencemethod）,也就是用差商代替微分方程及边界条件中的导数,最终化为代数方程求解；
（2）有限元法（finiteelementmethod）；
（3）把边值问題转化为初值问题，然后用求初值问题的方法求解。
2差分法
2.1—类特殊类型二阶线性常微分方程的边值问题的差分法
设二阶线性常微分方程的边值问题为
（&
2」）
（8.2.2）
yXx）-q（x）y（x）=f（x\a<
b,y（a）=a.y（Z?
）=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=<
x{<
•••<
畑<
xN=b
其中分点兀・=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<
/0<
^-1）,使y.=Mt而儿日和儿中至少有一个是小于M的.此时
因为録no,Mno,所以/（>
o,这与假设矛盾，故M只能是儿或片・证毕!
推论差分方程组（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<
b-
显然当/?
T0时，£
=y（兀•）—”t0.这表明当"
tO时，差分方程组（&
2.7）或
（8.2.8）的解收敛到原边值问题（8.7.1）,（8.7.2）的解.
例&
2.1取步长力=0.1,用差分法解边值问題
y*_4y=3x,0<
x<
1,
\，（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,
<
X_1-（2+4x0.12）yf+yf+l=3x0.12x/,i=2,3,…,&
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
兀
准确值y（x.）
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
从表8.2・1可以看出.差分方法的计算结果的精度还是比较高的.若要得到更精确的数值解，可用缩小步长力的方法来实现.
8.2.2一般二阶线性常微分方程边值问题的差分法
对一般的二阶微分方程边值问题
y\x）+p（x）yXx）+q（x）y（x）=f（x\a<
b,
*+a2y\a）-a,（8.2.12）
0ye）+02）'
O）=0，
假定其解存在唯一.
为求解的近似值，类似于前面的做法，
（i）把区间/等分，即得到区间I=[a,b]的一个网格剖分：
a=x0<
x,<
•••<
%、<
xn=b,
其中分点兀•=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,用差分法求下列边值问题的近似解，并将结果与精确解进行比较.
（x）=y\x）+2y（x）+cosx,0<
^/2,*y（0）=-0.3,y（/r/2）=-0.1,
精确解是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
准确值y（兀）
-0.3
-0.3
/T/16
-0.3137967
-0.3137446
2^/16
-0.3154982
-0.3154322
3^/16
-0.3050494
-0.3049979
4^/16
-0・2828621
-0・2828427
5龙/16
-0.2497999
-0.2498180
6^/16
-0・2071465
-0・2071930
7^/16
-0.1565577
-0.1566056
71/2
-0・1000000
-0.1000000
8.3有限元法
有限元法（finiteelementmethod）是求解微分方程定解问题的有效方法之一，它特别适用在几何、物理上比较复杂的问题.有限元法首先成功地应用于结构力学和固体力学，以后又应用于流体力学、物理学和其他工程科学.为简明起见，本节以线性两点边值问题为例介绍有限元法.
考虑线性两点边值问題
厶y=-（p（x）y'
（x））'
+g（x）y（x）=f（x）,a<
by（8.3.1）
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<
h0.
显然/]（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即
心）-/（刃>
只有当“（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<
xx<
-<
x{v・-<
xn=b.
区间厶=［和兀］称为单元,
长度勺=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）代入
升级会员 