常微分方程边值问题的数值解法.docx
- 文档编号:10062948
- 上传时间:2023-02-08
- 格式:DOCX
- 页数:28
- 大小:113.86KB
常微分方程边值问题的数值解法.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),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)代入
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 微分方程 边值问题 数值 解法