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

第8章常微分方程边值问题的数值解法81引言第7章介绍了求解常微分方程初值问题的常用的数值方法；本章将介绍常微分方程的边值问题的数值方法。只含边界条件BOUNDARYVALUECONDITION作为定解条件的常微分方程求解问题称为常微分方程的边值问题BOUNDARYVALUEPROBLEM为简明起见，我们以二阶边值问题为例介绍常用的数值方法。一般的二阶常微分方程边值问题BOUNDARYVALUEPROBLEMSFORSECONDORDERORDINARYDIFFERENTIALEQUATIONS为,811,YXFYXAXB其边界条件为下列三种情况之一1第一类边界条件THEFIRSTTYPEBOUNDARYCONDITIONS,YB2第二类边界条件THESECONDTYPEBOUNDARYCONDITIONSAY3第三类边界条件THETHIRDTYPEBOUNDARYCONDITIONS0101,YAYB000,定理811设811中的函数及其偏导数,在,FXYYFX,YFX,DXA上连续若1对所有，有；,Y,0YFX2存在常数，对所有，有,MD,YFXM则边值问题811有唯一解。推论若线性边值问题812,YXPYXQFXABAB满足1和在上连续；,PXQFX,2在上,B0则边值问题811有唯一解。求边值问题的近似解，有三类基本方法1差分法DIFFERENCEMETHOD，也就是用差商代替微分方程及边界条件中的导数，最终化为代数方程求解；2有限元法FINITEELEMENTMETHOD；3把边值问题转化为初值问题，然后用求初值问题的方法求解。82差分法821一类特殊类型二阶线性常微分方程的边值问题的差分法设二阶线性常微分方程的边值问题为821,YXQYFXABAB其中在上连续，且,QXF,AB0用差分法解微分方程边值问题的过程是I把求解区间分成若干个等距或不等距的小区间，称之为单元；,II构造逼近微分方程边值问题的差分格式构造差分格式的方法有差分法,积分插值法及变分插值法；本节采用差分法构造差分格式；III讨论差分解存在的唯一性、收敛性及稳定性；最后求解差分方程现在来建立相应于二阶线性常微分方程的边值问题821,822的差分方程I把区间等分，即得到区间的一个网格剖分,IABN,IAB,011NXX其中分点，并称之为网格节点GRIDNODES；步长,IXHIBANHII将二阶常微分方程822在节点处离散化在内部节点IX处用数值微分公式1,2IXN823241112,IIIIIIIIYXYXHYXH代替方程822中，得I,824112IIIIIIIYXYXQYXFRXH其中24IIHR当充分小时，略去式824中的，便得到方程821的近似方程IRX,825112IIIIIYQYFH其中,分别是的近似值,称式,IIIIQXF11,II11IIIXYX825为差分方程DIFFERENCEEQUATION，而称为差分方程825逼近方程822的截IR断误差TRUNCATIONERROR边界条件872写成8260,NY于是方程825,826合在一起就是关于个未知量，以及个方101,NY1程式的线性方程组827221122211,IIIIINNNQHYFHF这个方程组就称为逼近边值问题821,822的差分方程组SYSTEMOFDIFFERENCEEQUATIONS或差分格式DIFFERENCESCHEME，写成矩阵形式2211122233322211111NNNYQHHFQHFYQHHFQ828用第2章介绍的解三对角方程组的追赶法求解差分方程组827或828,其解称为边值问题821,822的差分解DIFFERENCESOLUTION由于825是用二01,NY阶中心差商代替方程821中的二阶微商得到的，所以也称式827为中心差分格式CENTEREDDIFFERENCESCHEMEIII讨论差分方程组827或828的解是否收敛到边值问题821,822的解，估计误差对于差分方程组827，我们自然关心它是否有唯一解；此外，当网格无限加密，或当时，差分解是否收敛到微分方程的解为此介绍下列极值原理0HIYIYX定理821极值原理设是给定的一组不全相等的数，设01,N8292,0,1,2IIIIIILQNH1若,则中非负的最大值只能是或；,IY0NIY0Y2若,则中非正的最小值只能是或01ILIN证只证1的情形，而2的情形可类似证明IYILY用反证法记，假设,且在中达到因为不全0MAXIINM0121,NYIY相等，所以总可以找到某个，使，而和中至少有一个01IN0IYM01IY0I是小于的此时M000000112IIIIIIIYLQH因为，所以,这与假设矛盾，故只能是或证毕0,IQ0ILYM0YN推论差分方程组827或828的解存在且唯一证明只要证明齐次方程组821011200,1,2,IIIIIINYLQYH只有零解就可以了由定理871知，上述齐次方程组的解的非负的最大值01,NY和非正的最小值只能是或而，于是证毕0Y0,NY,2I利用定理821还可以证明差分解的收敛性及误差估计这里只给出结果定理822设是差分方程组827的解，而是边值问题821,822的解IIYX在上的值，其中则有YXI,18211224,96IIIMHYXBA其中44MAXBMY显然当时，这表明当时，差分方程组827或0H0IIIXY0H828的解收敛到原边值问题871,872的解例821取步长，用差分法解边值问题143,01,0YX并将结果与精确解进行比较2234XYEE解因为，,由式827得差分格式1NH4,3QFX22111228900,2,38,439IIIIYY,其结果列于表82101Y,IXHI表821IIXIY准确值IYX0100101003329230033365620200649163006506043030093136900933461404011608310116348250501316725013197966060137528801378578707013088630131208780801084793010875539090066411400665865101000从表821可以看出,差分方法的计算结果的精度还是比较高的若要得到更精确的数值解，可用缩小步长的方法来实现H822一般二阶线性常微分方程边值问题的差分法对一般的二阶微分方程边值问题821212,YXPYXQFXABAB假定其解存在唯一为求解的近似值，类似于前面的做法，I把区间等分，即得到区间的一个网格剖分,IAN,IAB,011NXX其中分点，步长,IXHIBAHII对式8212中的二阶导数仍用数值微分公式241112,IIIIIIIIYXYXYXH代替，而对一阶导数，为了保证略去的逼近误差为，则用3点数值微分公式；另外2OH为了保证内插，在2个端点所用的3点数值微分公式与内网格点所用的公式不同，即8213211120000221,1,634,3,IIIIIIINNNNNYXHYXXYXYXHH略去误差，并用的近似值代替，便得III,IIIIIIPXQFX到差分方程组82141112002211,1,2,34,3,IIIIIIINNYYYFNHHYYH其中,是的近似值整理得,1IIIIIIQXPXFNIYIX8215120212212212134,1,IIIINNNHYQYPFN解差分方程组8215，便得边值问题8212的差分解01,Y特别地,若，则式8212中的边界条件是第一类边值条1212,0,件此时方程组7716为,YAB8216211211212111,32,IIIINNNNHQYHPYFHPPYFQ方程组8216是三对角方程组，用第2章介绍的解三对角方程组的追赶法求解差分方程组8216，便得边值问题8212的差分解01,YIII讨论差分方程组8216的解是否收敛到微分方程的解，估计误差这里就不再详细介绍例822取步长，用差分法求下列边值问题的近似解，并将结果与精确解/16H进行比较2COS,02,03,1YXYXX精确解是1SINCOSYX解因为，,由式8217得28NH,2,COSPXQFX差分格式2122112212661COS1603,16S,2,6,261CO7IIINNYYYYY2610,其结果列于表822083,01Y0,9IXHI表822IIIY准确值IYX0003031/16031379670313744622/16031549820315432233/16030504940304997944/16028286210282842755/16024979990249818066/16020714650207193077/1601565577015660568/2010000000100000083有限元法有限元法FINITEELEMENTMETHOD是求解微分方程定解问题的有效方法之一，它特别适用在几何、物理上比较复杂的问题有限元法首先成功地应用于结构力学和固体力学，以后又应用于流体力学、物理学和其他工程科学为简明起见，本节以线性两点边值问题为例介绍有限元法考虑线性两点边值问题8312,LYPXQXYFXABAB其中,10,CPXQ,FB此微分方程描述了长度为的可变交叉截面表示为的横梁在应力和QXPX下的偏差FY831等价性定理记,引进积分221C,C,ABXABYB83322DBAIYPXYQXYFXY任取，就有一个积分值与之对应，因此是一个泛函21C,YXIIFUNCTIONAL，即函数的函数因为这里是的二次函数，因此称为二次泛函,对泛函833有如下变分问题VARIATIONPROBLEM求函数，使得对任意21C,YAB,均有21YAB,834IY即在处达到极小,并称为变分问题834的解I可以证明定理831等价性定理是边值问题831,832的解的充分必要条件是使泛函YY在上达到极小，即是变分问题834在上的解IY21C,AB21C,AB证充分性设是变分问题的解；即使泛函在上21C,ABIYYI21C,AB达到极小，证明必是边值问题831,832的解Y设是任意一个满足的函数，则函数X2,0，21,XYXAB其中为参数因为使得达到极小，所以，即积分YIIYIY22BAIYPQFXXD作为的函数，在处取极小值，故0I8350DY计算上式，得02208DD2DBABAIYPXXQYXFXYXYPXQXYFX36利用分部积分法计算积分DDD,BBAABABAPYXXPYX代入式836，得0837D2D0BAIYPXYQXYFX因为是任意函数，所以必有X8380F否则，若在上某点处有,AB0X，000PYQXYFX不妨设，000XF则由函数的连续性知，在包含的某一区间上有,ABPYQXYF作02200,ABXAXBX显然，且，但2C,XABDPXYQXYFX，00BAF这与式837矛盾于是式838成立，即变分问题834的解满足微分方程831,且Y故它是边值问题831,832的解,YA必要性设是边值问题831,832的解，证明是变分问题834的解；YX即使泛函在上达到极小I21C,AB因为满足方程831，所以0PXYQXYF设任意，则函数满足条件，21,YXAB且于是2,X22222DD2DBABAAIIIYPYXQXFXYXFYXYPQX22DBBAPYQXFXXX因为，所以当时，,0,0X22D0BAPXQX即0IY只有当时，这说明使泛函在上达到极小证毕0X0IYIY21C,AB定理832边值问题831,832存在唯一解证明用反证法若都是边值问题831,832的解，且不相等，则由12,YX定理831知，它们都使泛函在上达到极小，因而I1C,AB且，221IY矛盾因此边值问题831,832的解是唯一的由边值问题解的唯一性，不难推出边值问题831,832解的存在性留给读者自行推导832有限元法等价性定理说明，边值问题831,832的解可化为变分问题834的求解问题有限元法就是求变分问题近似解的一种有效方法下面给出其解题过程第1步对求解区间进行网格剖分01,INAXXB区间称为单元，长度称为步长，若1,IIIIX,2IIH1MAXINH，则称上述网格剖分为均匀剖分,2IHN给定剖分后，泛函833可以写成22DBAIYPXYQXYFXY8391221IINNIXFXS记第2步构造试探函数空间。为了计算积分839，最简单的近似方法是将分段线性函数的集合作为试探函数空间。设分别是边值问题831,832的解在01,NYYX节点处的值，用分段线性插值01,NX8310111,IIIIIIXXYYIXH近似，并称式8310为单元形状函数ELEMENTSHAPEFUNCTION1,IIIYXI，为了将线性插值8310标准化，令，1IXTH则将变到轴上的单元记，则式8310可写成1,IIIIXT0,01,NTT831101,0,1IIYNTTY第3步建立有限元方程组将式8310代入泛函839，有1122DDIIIINNXXIYPQXFYX由式8311知212110101010121112DNIIIIIIIINIIIIIIIIIIIIIIIIIYIYXTHHXTNTYTTFTNTTHPXHQXYT12DIIIIIIIITT83121102NIIIIIIHFXTTYT式8312右端第1个求和号内的第项对应第个单元是关于的二次型，可以写成1,IIY,8313TIIYK其中，,8314T1,IIIY1,1,IIIII称为单元刚度矩阵ELEMENTSTIFFNESSMATRIX，IK831512,110,11,110D,IIIIIIIIIIIIIIIIIIIHPXTHQXTTKTTHT而式8312的第2个求和号内的第项可以写成8316TIIBY其中，T12,IIB1,IIIYY1200DDIIIIIIIHFXTTHFXT于是式8312中求和号内的项可以写成8317TTIIIIISKYB再令及矩阵T01,NYY2100,1III列列C则IIY于是式8317又可写成TT2IIIIIISCYKBCY8318IIII把式8318代入式8312右端求和号内，得8319TTT1112NNNIIIIIIIIYSYCKYCBY记1,1,T11000000IINNIIIIIIKKC，11002213312442311,2,1,1,NNNNKKKK8320TT1211,0,0,NNIIIIBBCB,83212341T22,NNB则式8319化为8322TT2IYYKB并称为总刚度矩阵TOTALSTIFFNESSMATRIX，称为右端向量K因为是使取极小值的函数，所求的自然使式8322的右边取YXI01,NY极小值，所以应有8323TTD2,2,IIYYKB记,则式8823为12,IJNNKB0000D2NRJRRRJINNRIIJIJYBYYK02,1,2,NIRIYBIN或83240,0,NIRIKY得方程组8325YB因为是已知的，不能任意选取，所以不能要求式8323对也成立0,NY0,NY因此方程组8324或8325中应当去掉首末2个方程，并把其他方程中含有的改为已知量，所得方程组就是未知量满足的代数方程组11,NY122113322211,3,2,2,1,NNNNNYKKKYK832612103211,NNBBK方程组8325或8326称为有限元方程组FINITEELEMENTSYSTEMOFEQUATIONS用第2章介绍的解三对角方程组的追赶法求解有限元方程组8326，可解出，即得变分问题834的近似解，也就是边值问题831,832解的近似值11,NY这种方法称为有限单元法FINITEELEMENTMETHOD或有限元法例831用有限元法求下列边值问题的近似解，并将结果与精确解进行比较2481,0,01,XYXX取，精确解是02H2解因为，,由式8326得有限5NH2,4PXQFX元方程组12343670867,4363567,852YYYY其结果列于表831表831IIXIY准确值IYX0000102016440162040244302430602434024408016200165100上面虽然是对边值问题831,832介绍的有限元解法，但其解题步骤对一般的微分方程定解问题也是适用的对所给微分方程定解问题，首先找出相应微分方程的变分形式，然后进行下列步骤第1步对定义区域或定义区间进行网格剖分,将定义区域或定义区间剖分为若干个小单元一维是区间；多维是区域，如矩形、三角形等；第2步构造试探空间即选择适当的插值函数类第3步建立有限元方程组计算单元刚度矩阵及右端向量，再形成总刚度矩阵及总右端项，写出有限元方程组；结合定解条件，求解有限元方程组注从形式上看，用有限元法解微分方程定解问题很繁，但是有限元法的求解步骤规范，便于在计算机上实现，并且总刚度矩阵是三对角对称正定矩阵，因此有限元方程组可用第2章介绍的解三对角方程组的追赶法求解有限元法最主要的优点是可以处理相当复杂的区域上的边值问题。84打靶法解常微分方程边值问题的打靶法SHOOTINGMETHOD，也称为尝试法，其基本思想是把边值问题转化为初值问题来解首先作出一些只满足一端边值条件的解，然后再从这些解中找出适合另一端边值条件的解下面以二阶常微分方程带第一类边界条件的边值问题8412,YXFYXAXBAB为例介绍常微分方程边值问题的打靶法76节曾介绍过二阶常微分方程初值问题7610,YTFTYATBA的求解方法将上式中的改为，改为，得TXS843,YFYXX设初值问题843的解为，显然依赖于，即为了求解边值1X1S1,YS问题841,842，必须求，使之满足S,YB下面介绍2种方法来求方法1根据实际问题情况选一个，解初值问题1S8441,YFXYAXA得到的解仍记为1X若或为事先给定的精度，则就是边值问题841,1YBY1YX842的解否则，根据与的误差，将修改为；例如取,再解初值问1B1S221S题8452,YXFYXAXBAS得到解2YX若满足或，则就是边值问题841,842的解B222YX否则，再将适当修改为例如可用作线性插值求S3S1,S3S,211然后解初值问题8463,YXFYXAXBAS的解如此继续下去，直到达到精度要求为止方法2求另一种方法是求函数的一个零点S,FYBS对于每一个自变量，通过解初值问题844，可解出计算,YX,847,S然后采用第3章介绍的求方程根的方法求的零点首先寻找，使，则在区间或12,S120,F12,S内至少有的一个零点然后可用二分法求也可用NEWTON迭代公式2,S1KKFSS求的近似值S几何解释微分方程边值问题841,842的解是一条通过YX两点的曲线如图841初值问题844的解是一条通过,AAYBBY1,S点、斜率为的曲线如图841初值问题845的解是一条1AS2通过点、斜率为的曲线如图841这有点象射击者从定点向目,2A标射击一样根据经验以某一角度斜率为试射一次如果与目标相差太大，调整射1S击角度斜率为，再进行射击如此继续进行下去，直到击中目标，或击中的点与的2SB误差在允许的范围之内。图841参考文献1提供的资料还讨论了选择初始值的重要性，并介绍了多重打靶YAS法，这里就不作详细介绍，有兴趣的读者可参看相关资料85算法程序851用差分法解二阶线性常微分方程的边值问题用差分法解二阶线性常微分方程的边值问题A,B是区间,STEP是步长,ALPHA,BETA是初值F,Q分别是式821中的FX,QXFUNCTIONYDIFFMETHODF,Q,A,B,STEP,ALPHA,BETANBA/STEPXASTEPBAZEROSN1FORI1N1AI,I12FEVALQ,XI1STEP2IFIN1AI,I11AI1,I1ENDENDB1STEP2FEVALF,X2ALPHAB2N2STEP2FEVALF,X3N1BN1STEP2FEVALF,XNBETABBL,ULUAYULBFORI1LENGTHYSPRINTF107F,YIENDEND例851取，用差分法求下列边值问题的近似解01H4,01,0,12YXX解在MATLAB命令窗口输入FINLINE4XQINLINE4DIFFMETHODF,Q,0,1,01,0,2回车，可的结果。852用差分法解一般二阶常微分方程的边值问题用差分法解一般二阶常微分方程的边值问题A,B是区间,STEP是步长,ALPHA,BETA是初值F,P,Q分别是式8212中的FX,PX,QXFUNCTIONYDIFFMETHOD_2F,P,Q,A,B,STEP,ALPHA,BETANBA/STEPXASTEPBAZEROSN1A1,122STEP2FEVALP,X2A1,22STEPFEVALF,X3FORI2N2AI,I22FEVALP,XI1STEP2AI,I12STEPFEVALF,XI1AI1,I2STEPFEVALF,XI1ENDAN1,N22STEPFEVALF,XNAN1,N122FEVALP,XNSTEP2B12STEP2FEVALQ,X22STEPFEVALF,X2ALPHAB2N22STEP2FEVALQ,X3N1BN12STEP2FEVALQ,XN2STEPFEVALF,XNBETABBL,ULUAYULBFORI1LENGTHYSPRINTF107F,YIENDEND例852取，用差分法求下列边值问题的近似解/16H2COS,02,03,1YXYXX解在MATLAB命令窗口输入FINLINECOSXPINLINE1QINLINE2DIFFMETHOD_2F,P,Q,0,PI/2,PI/16,03,01回车，可得结果。853用有限元法解二阶常微分方程的边值问题用有限元法解二阶常微分方程的边值问题A,B是区间,H是步长,ALPHA,BETA是初值F,P,Q分别是式831中的FX,PX,QXFUNCTIONYFINIELEMF,P,Q,A,B,STEP,ALPHA,BETANLENGTHSTEP1X1AFORI1N1XI1XISTEPIENDSYMSTFF1/STEPN1FEVALF,XNTSTEPN1STEPN1FEVALP,XNTSTEPN11TTKNNNINTFF,T,0,1SYMSTFF1/STEP2FEVALF,X1TSTEP2STEP2FEVALP,X1TSTEP21TTK101INTFF,T,0,