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

1、第8章 常微分方程边值问题的数值解法第7章介绍了求解常微分方程初值问题的常用的数值方法；本章将介绍常微分方程的边值问题的数值方法。只含边界条件(boundary-value condition)作为定解条件的常微分方程求解问题称为常微 分方程的边值问题(boundary-value problem).为简明起见，我们以二阶边值问题为例介绍常 用的数值方法。一般的二阶常微分方程边值问题(boundary-value problems for second-order ordinarydifferential equations)为y (x) f (x, y(x), y (x), a x b,(8

2、.1.1)其边界条件为下列三种情况之一：(1)第一类边界条件(the first-type boundary conditions) ： y(a) , y(b) ;(2)第二类边界条件(the second-type boundary conditions) ： y (a), y (b);(3)第三类边界条件(the third-type boundary conditions):y(a)0y (a)1，y(b) 0y (b)1，° 0,° 0,000.定理8.1.1 设中的函数f (x, y, y )及其偏导数fy(x, y,y), fy(x, y, y )在D (x,

3、y,y) a x b, y , y 上连续.若(1)对所有(x, y, y) D,有 fy(x,y,y) 0；(2)存在常数 M ,对所有(x, y, y) D ,有 fy (x, y, y ) M ,则边值问题(8.1.1)有唯一解。推论若线性边值问题(8.1.2)y (x) p(x)y (x) q(x)y(x) f (x), a x b, y(a) , y(b)满足(1) p(x), q(x)和 f (x)在a, b上连续;(2)在a,b上，q(x) 0,则边值问题(8.1.1)有唯一解。求边值问题的近似解，有三类基本方法：(1)差分法(difference method),也就是用差商代

4、替微分方程及边界条件中的导数，最终化为代数方程求解；(2)有限元法(finite element method);(3)把边值问题转化为初值问题，然后用求初值问题的方法求解。差分法8.2.1一类特殊类型二阶线性常微分方程的边值问题的差分法设二阶线性常微分方程的边值问题为y (x) q(x)y(x) f(x), a x b,(8.2.1)y(a) , y(b) ,(8.2.2)其中q(x), f(x)在a,b上连续,且q(x) 0.用差分法解微分方程边值问题的过程是：(i)把求解区间a,b分成若干个等距或不等距的小区间，称之为单元；(ii)构造逼近微分方程边值问题的差分格式.构造差分格式的方法有

5、差分法，积分插值法及变分插值法；本节采用差分法构造差分格式；(iii)讨论差分解存在的唯一性、收敛性及稳定性；最后求解差分方程.现在来建立相应于二阶线性常微分方程的边值问题(8.2.1)，的差分方程.把区间I a,b N等分，即得到区间I a,b的一个网格剖分: a x0 Xi LXn 1 Xn b ,其中分点xi a ih (i 0,1,L , N ),并称之为 网格节点(grid nodes);步长h(ii )将二阶常微分方程(8.2.2)在节点为处离散化：在内部节点Xi (i 1,2,L ,N 1)处用数值微分公式y (Xi)y(K 1) 2y(x) y(Xi)h1。()12Xi 1i

6、Xi(8.2.3)代替方程(8.2.2)中y (xi),得(8.2.4)y(" 1)2y(2X) y(" 1)q(Xi)y(x) f(x) R(x), h2其中"x) -y()(i).当h充分小时，略去式(8.2.4)中的R(x),便得到方程的近似方程yi 1 2yiyi 1h2qiyifi,(8.2.5)其中 qi q(Xi), fif (x), x 1, yi, yi 1 分别是 y(x 1), y(x) y(x 1)的近似值 , 称式(8.2.5)为差分方程(difference equation),而R(x)称为差分方程逼近方程的截断误差(truncati

7、on error).边界条件写成y。Yn(8.2.6)于是方程(8.2.5),合在一起就是关于N1个未知量y0,y1,L , yN，以及N1个方程式的线性方程组： 22 -(2 qih )yiy2h fl,2、. 2y i(2 qh)y yihfi, i 1,2,L , N 1,(8.2.7)yN 2(2 qN 1h )yN1hfN1这个方程组就称为逼近边值问题(8.2.1),的差分方程组(system of difference equations)或差分格式(difference scheme),写成矩阵形式(2 q1h2)11(2 q2h2)11(2 q3h2)1O O1(2y1h2 f

8、1y2h2f2y3h2f3OMM,22 rqN 2h )1yN 2 h fN 2221(2 qN 1h ) yN 1 h fN 1(8.2.8)用第2章介绍的解三对角方程组的追赶法求解差分方程组(8.2.7)或 其解y0,y1,L , yN称为边值问题的差分解(difference solution).由于是用二阶中心差商代替方程中的二阶微商得 到的，所以也称式为 中心差分格式(centered-difference scheme).(iii )讨论差分方程组(8.2.7)或的解是否收敛到边值问题的解，估计误差.对于差分方程组(8.2.7),我们自然关心它是否有唯一解；此外，当网格无限加密，或

9、当h 0时，差分解yi是否收敛到微分方程的解y(x”.为此介绍下列极值原理：定理8.2.1 (极值原理)设y0, y1,L "n是给定的一组不全相等的数，设l( )y-2yiyqqi 0, i 1,2,L,N.(8.2.9)hN(1)若l(yi)0, i1,2,L , N ,则j0中非负的最大值只能是y0或yN ;N , 一一一(2)若l(yi)0, i1,2,L ,N ,则yii0中非正的最小值只能是y0或yN.证 只证(1) l(yj 0的情形，而(2) l(yi)0的情形可类似证明.i0 (1 io N 1),使yi0M ,而yi。1和yi01中至少有一个用反证法.记M max

10、 yi ，假设M 0 ,且在y1, y2,L , yN 1中达到.因为yi不全 相等，所以总可以找到某个 是小于M的.此时i(yjyio i2%。h2yio iqio yioM 2Mh2qi°MqiM因为q0 0, M 0,所以l(yi。) 0,这与假设矛盾，故 M只能是y。或yN.证毕！ 推论 差分方程组(8.2.7)或的解存在且唯一.证明 只要证明齐次方程组l(Yi) yi1 2yi yi1 qiYi 0, q 0, i 1,2,L ,N 1,h(8.2.10)y0 0, yN 0只有零解就可以了 .由定理8.7.1知，上述齐次方程组的解y0,y1,L ,yN的非负的最大值和非正

11、的最小值只能是y0或yN.而y00, yN 0 ,于是y 0, i 1,2,L ,N.证毕！利用定理8.2.1还可以证明差分解的收敛性及误差估计.这里只给出结果：定理8.2.2设yi是差分方程组的解，而 y(x)是边值问题 的解y(x)在Xi上的值，其 中i 0,1,L ,N .则有M4h22iy(Xi) y (b a) ,(8.2.11)96其中 M 4 max y(x). a x b显然当h 0时，i y(x) yi 0.这表明当h 0时，差分方程组(8.2.7)或的解 收敛到原边值问题的解.例8.2.1 取步长h 0.1 ,用差分法解边值问题y 4y 3x,0x1,y(0) y(1) 0

12、,并将结果与精确解 y(x) 3 e2x e2x/4e2 e23x4进行比较.解 因为N1/h 10, q(x) 4, f(x)3x,由式(8.2.7)得差分格式(2 4 0.12)y1y2 3 0.12 0.1,y 1 (2 4 0.1% y 1 3 0.12x, i 2,3,L ,8,_2_2_y (2 4 0.12)y9 3 0.12 0.9,y0y100, x 0 ih 0.1i, i 1,2,L ,9 ,其结果列于表 8.2.1.表 8.2.1ixiyi准确值y(xi)010010.0332923一2-0. 0649163一3-0. 0931369一4-0. 1160831一5-0.

13、 1316725一6-0. 1375288一7-0. 1308863一8-0. 10847931一9-0. 0664114一1000从表8.2.1可以看出，差分方法的计算结果的精度还是比较高的.若要得到更精确的数值解，可用缩小步长 h的方法来实现.8.2.2一般二阶线性常微分方程边值问题的差分法对一般的二阶微分方程边值问题y (x) p(x)y(x) q(x)y(x) f(x), a x b,iy(a) 2y (a),(8.2.12)iy(b) 2y(b), 假定其解存在唯一. 为求解的近似值，类似于前面的做法，(i )把区间I a,b N等分，即得到区间Ia,b的一个网格剖分：axo xi

14、Lxn i xn b ,其中分点 xia ih (i 0,1,L,N),步长 h bNa.(ii )对式(8.2.12)中的二阶导数仍用数值微分公式y (xi)y(xi 1) 2y(xj y(x 1)h"04)(i)，x12代替，而对一阶导数，为了保证略去的逼近误差为O(h2),则用3点数值微分公式；另外为了保证内插，在2个端点所用的3点数值微分公式与内网格点所用的公式不同，即2y(xi)i12h i1 5y ()x1 i xi, i 1,2，L ,N 11y (xo)23y(%) 4y(x) y(x?) h1y(o)，2hxo(8.2.13)y (xn)y(xN 2) 4y(xn

15、1) 3y(xn) h2h1y(N ),xNf(xi),便得略去误差，并用 y(xj的近似值yi代替y(xi), Pi p(x) qi q(x),片到差分方程组h2(yi12ylyio2h(yi1yi1)qiyifi,i 1,2,L,N1,iyo2( 3yo 4yi 、2,(8.2.14)2h1 yN(yN 2 4yN 1 3yN),2h其中qiq(Xi), pi p( xi), f if(xi),i 1,2,L,N 1, yi 是 y (xi)的似值.整理得(2h 1 3 2)yo 4 2Y12y22h ,2 2 _.(2 hpi)yi 1 2(2 h qi)yi (2 hpi)yi 1 2

16、h fi, i 1,2,L,N 1,(8.2.15)2 yN 24 2yN 1(3 2 2h 1)yN 2h解差分方程组(8215)，便得边值问题的差分解y0, y1,L , yN .特别地，若1 1, 20, 1 1, 20，则式(8.2.12)中的边界条件是第一类边值条件：y(a) , y(b) ;此时方程组为2(2 hi)% (2 hp'y22h2R (2 hpj ,(2 hpi)yi 1 2(2 h2qi)yi (2 hpi)yi 1 2h2fi, i 2,3,L,N 2,(8.2.16)22(2 hpN 1 )yN 2 2(2 h qN 1)yN 1 2h fN 1 (2 h

17、pN 1) .方程组(8.2.16)是三对角方程组，用第2章介绍的解三对角方程组的追赶法求解差分方程组，便得边值问题的差分解 y0, y1,L ,yN.(iii )讨论差分方程组(8.2.16)的解是否收敛到微分方程的解，估计误差.这里就不再详细介绍.例8.2.2 取步长h /16，用差分法求下列边值问题的近似解，并将结果与精确解进 行比较.y (x) y (x) 2y(x) cosx, 0 x 2, y(0)0.3, y( 2)0.1,1精确斛y(x) (sin x 3cos x).p(x) 1, q(x) 2,. 16( 1) Y2.16 2 cos 162.16 2 ( 2) yi2.1

18、6 2 cos i .16 ,16 2 ( 2) yN 116 2 cos 7 .162f (x) cosx,由式(8.2.17)得16 ( 1)( 0.3),16( 1) yi 1i 2,3,L ,6,.16 ( 1)( 0.1),10解因为 N ( /2 0)/h 8，差分格式2 2，/16 2 ( 2) y1222.16 ( 1) yi1 2 222.16 ( 1) Yn 2 2 22y00.3, y80.1,为 0 ih 0.1i, i 1,2,L ,9,其结果列于表 8.2.2.表 8.2.2ixiyi准确值y(xi)001/1622 /1633 /1644 /1655 /1666

19、/1677 /168/2有限元法有限元法(finite element method)是求解微分方程定解问题的有效方法之一，它特别适用在几何、物理上比较复杂的问题.有限元法首先成功地应用于结构力学和固体力学，以后又应用于流体力学、物理学和其他工程科学.为简明起见，本节以线性两点边值问题为例介绍有限元法.考虑线性两点边值问题Lyp(x)y(x)q(x)y(x) f (x), a x b,(8.3.1)y(a),y(b),(8.3.2)其中 p(x) 0, q(x) 0, p C1a,b, q, f Ca,b.此微分方程描述了长度为b a的可变交叉截面(表示为q(x)的横梁在应力p(x)和f(x)

20、下的偏差y(x).8.3.1 等价性定理记C2a,by y y(x)C2a,b, y(a) , y(b) ，引进积分,，、b2, 、2,I(y) a p(x)y(x) q(x)y (x) 2f (x)y(x) dx.(8.3.3)a2任取y y(x) C1a,b,就有一个积分值I(y)与之对应，因此I(y)是一个泛函(functional),即函数的函数.因为这里是y , y的二次函数，因此称 I(y)为二次泛函.对泛函(8.3.3)有如下变分问题(variation problem)：求函数 % C2a, b,使得对任意y C2a,b,均有I(y) I (%,(8.3.4)即I(y)在处达到

21、极小，并称为变分问题(8.3.4)的解.可以证明：定理8.3.1(等价性定理)是边值问题的解的充分必要条件是%使泛函I (y)在22到极小, 设其中2_ 一 一(充分性)设为 Cia,b是变分问题 证明为、是边值问题(8.3.1),的解.2(X)是C a,b任意一个满足(a)y(x) %x)I (y)的解；即y%泛函(b) 0的函数，则函数_2(x) Cia,b,为参数.因为%使彳#I(y)达到极小，所以1(%) 1(%,2.I(y)在 Cia,b上达即积分1(%22P(x)%(x)(x)q(x) %x)(x)2 f (x)%x)(x) dx作为的函数，在0处取极小值I (%,故计算上式，得F

22、i(%)1(%(8.3.5)bP(x)%(x) a(x)2q(x)%x)(x)2 2f(x)%x)(x)dx2 P(x)%(x)(x) (x)2q(x)yfx)(x) (x) 2f(x)(x) dxp(x)%(x) (x)q(x)%x)(x) f (x) (x) dx.(8.3.6)利用分部积分法计算积分P(x)%(x)(x)dxba P(x)%(x)d (x)P(x)%(x) (x)abb(x) P(x)%(x)adx代入式(8.3.6),a (x) P(x)%(x) dx.I(%b2a 0P(x)%(x)q(x)%x) f(x) (x)dx 0.(8.3.7)因为(x)是任意函数，所以必有

23、P(x)%(x)否则，若在a,b上某点X0处有q(x)%x) f (x) 0.(8.3.8)p(Xo)%(X0)qd)%。) f (Xo) 0,C2a,b上达到极小，即 为是变分问题在C2a,b上的解.不妨设P(Xo)%(Xo)q(Xo)%X。)f (Xo) 0,则由函数的连续性知，在包含 X0的某一区间a0,b0上有p(x)%(x)q(x)%x) f (x) 0.(x)0,x a,bao，bo,(x a°)2(x bo)2, ao x b0.显然(x) C2a,b,且(a)(b) 0,但p(x)%(x)q(x)%x) f (x) (x)dxb0a02- -设任息y y(x) Cia

24、,b,则函数 (x)2且(x) C a,b.于是I(y) 1(% 1(% )1(%b22. p(x)%(x)(x)q(x)%x)(x)y(x) %x)满足条件(a)2 f (x) %x)(x) dx(b) 0,2p(x)%(x)2q(x) %x)2 f (x)%x) dxb2 a p(x)%(x) (x) q(x)%x) (x) f(x) (x) dx a22p(x) (x)q(x) (x) dxb2p(x)%(x)aq(x)%x) f (x) (x)dxp(x) (x)2q(x) (x)2 dx22p(x) (x)q(x) (x) dx.a因为p(x) 0, q(x) 0,所以当 (x)b9

25、90时， p(x) (x)2 q(x) (x)2 dx 0,即 aI(y) I(%0.p(x) %(x)q(x)%x) f (x) (x)dx 0,这与式(8.3.7)矛盾.于是式成立，即变分问题的解例满足微分方程且 炖a) , %b)故它是边值问题的解.(必要性)设% %x)是边值问题(8.3.1),的解，证明 是变分问题的解；即：使泛2函I(y)在CJa,b上达到极小.因为 % %x)满足方程(8.3.1),所以 p(x)%(x)q(x)%x) f(x) 0.只有当(x) 0时，I(y) I(%0.这说明%吏泛函I(y)在C2a,b上达到极小.证毕!定理8.3.2边值问题存在唯一解证明 用

26、反证法.若y(x), y2(x)都是边值问题(8.3.1),的解，且不相等，则由定理知, 2它们都使泛函I (y)在C1 a,b上达到极小，因而I(y1) I(y2)且 I(y2) My),矛盾！因此边值问题(8.3.1),的解是唯一的.由边值问题解的唯一性，不难推出边值问题(8.3.1),解的存在性(留给读者自行推导).8.3.2 有限元法等价性定理说明，边值问题(8.3.1),的解可化为变分问题的求解问题.有限元法就是求变分问题近似解的一种有效方法.下面给出其解题过程：第1步 对求解区间进行网格剖分a x0x1 L x Lxnb,区间L 41, x 称为单元，长度 hi xi x 1 (i

27、 1,2,L ,n)称为步长，h maxhi .若 1 i nh hi (i 1,2,L ,n),则称上述网格剖分为均匀剖分给定剖分后，泛函(8.3.3)可以写成b22.I(y) a P(x)y(x) q(x)y (x) 2f (x)y(x) dxanxii 1 x1记n(8.3.9)P(x)y(x)2 q(x)y2(x) 2f (x)y(x)dx=§.i 1第2步构造试探函数空间。为了计算积分(8.3.9),最简单的近似方法是将分段线性函数的集合作为试探函数空间。设y0, y1,L ,yn分别是边值问题 的解y(x)在节点 Xo,Xi,L , xn处的值，用分段线性插值_xi x

28、x xi 1yy(x)'一yi 1 -x, x h ”为.(8.3.io)hihi近似 y(x),x IiXi 1, x ,并称式(8.3.10)为单元形状函数(element shape function)为了将线性插值(8.3.10)标准化，令hi则将Ii Xi1, Xi变到 t轴上的单元0,1.记 N°(t) 1 t, N1(t) yN0(t)yN1(t)yi, t 0,1.建立有限元方程组.将式(8.3.10)代入泛函，有t ,则式(8.3.10)可写成(8.3.11)I(y)I(y)nXi22P(x)y(x)q(x)y(x) dxf (x)y(x)dx .由式(8.

29、3.11)知I(y) 1p(Xiithi)(yiyi1)hiq(Xithi)N°(t)yNi(t)yJ2dti i 0hin i20hf(xi thi)N0(t)yi i Ni(t)yjdti 1n i0 hi ip(Xi i thi) hiq(x i thi)(i t)2 y2i i i一.i ,.、.，2hip(Xii th)hq(xi th)(i t)t ihii p( i thi) hiq(X i thi)t2 y2dtn i20hif(x i thi)(i t)yi i tyi dt.(8.3.i2)i 1式(8.3.12)右端第1个求和号内的第i项(对应第i个单元)是关于

30、yi 1,yi的二次型，可以写成(8.3.13)(Y)TKY其中y(yi i,yi)T, K(i)Ki 1,i 1K0)Ki,i 1(i)Ki 1,iKi,i(8.3.14)K(i)称为单元刚度矩阵(element stiffness matrix), .i ._Ki 1,i i 0 hi p(x i thi)hq(K ithi)(1 t) dt,.(i)1 .12.Ki,i 0 hi p(Xi ith)hq(X ithi)tdt,(8.3.15)iKi 1 Kii i 0 hi p(Xi ith)hiq(%i -)(1 t)t dt.而式(8.3.12)的第2个求和号内的第i项可以写成(b)

31、TY(8.3.16)其中b(bi(i),b2i)T, Y(i)(yii,yi)T,(i)1 .(i)1 _bihi0 f(Xiithi)(1t)dt,b2hi0f(Xiithjtdt.于是式(8.3.12)中求和号内的项可以写成Si(Y )TK(i)Y 2(b)TY .(8.3.17)再令Y (y0,yi,L ,yn)T及2 (n 1)矩阵0L0100L00L0010L0则 Y(i) C(i)Y .于是式(8.3.17)又可写成Si(C (i)Y )TK (i)(C(i)Y) 2(b)T (C Y)(8.3.18)YT(C(i)TKC(i)Y2 (C(i)Tb(i) TY .把式(8.3.18

32、)代入式右端求和号内，得i(y)YTn(C(i)TKC Y(C)Tb(8.3.19)(C(i)TK(i)C(i)0(i)Ki 1,iKKi,i 100 (i)Ki 1,iKi(i)0M 0010K01(1) K11 K11K21K12(2) K22 K22八32K23 K33K33O(4) K34O(n 1)n 1,n 2(n 1) n 1,nO(n)1 n 1,n(n)n,n 1(n)n 1,n(n)Kn,n(8.3.20)b(C)Tbi 1n(0,Li 1i 1, i,0,b1 b2i)QL ,0)T (bi , b2bi,b22)b1(3) ,b23)bi4) ,L ,b2n 1)b1(

33、n),b2(n)T,(8.3.21)则式(8.3.19)化为I(y)YTKY 2bTY.(8.3.22)并称K为总刚度矩阵(total stiffness matrix),称b为右端向量.因为y(x)是彳I (y)取极小值的函数，所求的y0, y1,L , yn自然使式(8.3.22)的右边取极小值，所以应有d(8.3.23)一(YTKY 2bTY) 0, i 0,1,2,L,n. dyi记 K (Kj)nn,b(b1, b2,L ,bn)T ,则式(8.8.23)为d ndyi r 0nKrjYj j 0nyr2br yrr 0nKriYr r 0Ki jyj 02bn2 Kir yrr 0

34、0,1,2,L ,n 1,nKir Yr r 00,1,L ,n,(8.3.24)得方程组KY(8.3.25)是已知的，不能任意选取,所以不能要求式此方程组或中应当去掉首末 2个方程，并把其他方程中含有的(8.3.23)对yo, yn也成立.因 yo, yn改为已知量，所得方程组就是未知量y1,y2,L ,yn1满足的代数方程组 K11K11K21K12 K22K22OyiK23O(n 2)Kn 2,n 3(n 2)Kn 2,nO(n 1)2 Kn 2,n 2y2M(n 1)Kn 1,n 2(n 1)Kn 2,n 1(n 1)(n)Kn 1,n 1 Kn 1,nyn 2yn 1(1) (2)(

35、1)b2b1K10(2) (3)b2b1(8.3.26)M(n 2)(n 1)b2b1(n 1) (n)(n)b2b1Kn 1,n方程组(8.3.25)或称为 有限元方程组(finite element system of equations).用第2章介绍的解三对 角方程组的追赶法求解有限元方程组(8.3.26),可解出y1,y2,L , yn 1 ,即得变分问题的近似解，也就是边值问题解的近似值.这种方法称为 有限单元法(finite element method) 或有限元法.例8.3.1用有限元法求下列边值问题的近似解，并将结果与精确解进行比较2(xy(x) 4y(x) 4x 8x 1

36、, 0 x 1,y(0) y(i) 0,取h 0.2,精确解是y(x) x2 x.解因为 N 1/h 5, p(x) x, q(x) 4, f(x) 4x2 8x 1,由式(8.3.26)得有限元方程组2.53333y1 1.36667y20.082667,1.36667yl 4.53333y2 2.36667y30.306667,2.36667y2 6.53333y3 3.36667y40.466667,3.36667y3 8.53333y40.562667.其结果列于表8.3.1.表 8.3.1ixiyi准确值y(xi)0000123145100上面虽然是对边值问题(8.3.1),介绍的有

37、限元解法，但其解题步骤对一般的微分方程定 解问题也是适用的.对所给微分方程定解问题，首先找出相应微分方程的变分形式，然后进行下列步骤：第1步 对定义区域(或定义区间)进行网格剖分，将定义区域(或定义区间)剖分为若干 个小单元(一维是区间；多维是区域，如矩形、三角形等)；第2步构造试探空间；即选择适当的插值函数类.第3步 建立有限元方程组.计算单元刚度矩阵及右端向量，再形成总刚度矩阵及总右端项，写出有限元方程组；结合定解条件，求解有限元方程组注从形式上看，用有限元法解微分方程定解问题很繁，但是有限元法的求解步骤规范，便于在计算机上实现，并且总刚度矩阵是三对角对称正定矩阵，因此有限元方程组可用第2

38、章介绍的解三对角方程组的追赶法求解.有限元法最主要的优点是可以处理相当复杂的区域上的边值问题。打靶法解常微分方程边值问题的 打靶法(shooting method),也称为尝试法，其基本思想是把边 值问题转化为初值问题来解： 首先作出一些只满足一端边值条件的解，然后再从这些解中找出适合另一端边值条件的解.下面以二阶常微分方程带第一类边界条件的边值问题(8.4.7)(8.4.2)y(x) f (x,y(x),y (x), a x b y(a) , y(b).为例介绍常微分方程边值问题的打靶法.节曾介绍过二阶常微分方程初值问题(7.6.10)y (t) y(a)的求解方法.将上式中的t改为x ,y

39、 (x) y(a)f(t, y(t),y(t), a t b, ,y(a),改为s ,得f (x,y(x),y (x), a x b, ,y (a) s,(8.4.3)设初值问题(8.4.3)的解为y1 (x),显然y1(x)依赖于s ,即y1(x)y1(x, s).为了求解边值问题，必须求s % 使之满足y1 (b)y1 (b, %.下面介绍2种方法来求能方法1根据实际问题情况选一个s1,解初值问题(8.4.4)y f (x, y, y), a x b, y(a) , y(a) 5,得到的解仍记为y1(x).若yi(b)或y1(b)(为事先给定的精度)，则yi(x)就是边值问题(8.4.1)

40、,的解.否则，根据1y1(b)与 的误差，将s修改为电；例如取s 一s1，再解初值问题1y (x) f (x,y(x),y (x), a x b,(8.4.5)y(a) , y(a)s2,得到解y2(x).若2丫29)满足2 或2，则y2(x)就是边值问题(8.4.1),的解.否则，再将 s适当修改为s3;例如可用(1,s1), ( 2,s2)作线性插值求s3：s ss2 s1 ()s3 sl(1 ) ,21然后解初值问题y (x) f (x,y(x),y (x), a x b,(8.4.6)y(a) , y(a)s3,的解.如此继续下去，直到达到精度要求为止方法2求能另一种方法是求函数F(s

41、) y(b, s)的一个零点 能对于每一个自变量 s ,通过解初值问题(8.4.4),可解出y(x,s).计算F(s) y(b,s)然后采用第3章介绍的求方程根的方法求 F (S)的零点.首先寻找 s(), s(),使 F(s() 0, F(s ) 0 ,则在区间(s(), s()或(s(), s()内至少有F (s)的一个零点.然后可用二分法求 %也可用Newton迭代公式sk 1skF(Sk)F回)求哪近似值.几何解释：微分方程边值问题(8.4.1),的解y(x)是一条通过A(a, y(a), B(b, y(b)两 点的曲线(如图 初值问题的解y(x,s1)是一条通过点 A(a,y(a)、

42、斜率为y (a) s的曲线 (如图 初值问题的解y(x, S2)是一条通过点 A(a,y(a)、斜率为y (a) S2的曲线(如图 这 有点象射击者从定点 A向目标B射击一样.根据经验以某一角度(斜率为、)试射一次.如 果与目标相差太大，调整射击角度(斜率为s2),再进行射击.如此继续进行下去，直到击中 目标，或击中的点与 B的误差在允许的范围之内。并介绍了多重打靶法，这里就不作详细介绍，有兴趣的读者可参看相关资料算法程序8.5.1 用差分法解二阶线性常微分方程的边值问题%用差分法解二阶线性常微分方程的边值问题%a,b是区间,Step是步长,Alpha, Beta是初值% f, q 分别是式(

43、8.2.1)中的 f (x), q (x)function Y = DiffMethod(f, q, a, b, Step, Alpha, Beta)N = (b-a) / Step;X = a : Step : b;A = zeros( N-1 );for i = 1 : N-1A(i,i) = -1 * ( 2+feval(q, X(i+1) * StepA2 ); if i = N-1A(i,i+1) = 1;A(i+1,i) = 1; endendB(1) = StepA2 * feval(f, X(2) - Alpha;B(2:N-2) = StepA2 * feval(f, X(3

44、:N-1);B(N-1) = StepA2 * feval(f, X(N) - Beta;B = B'L,U = lu( A );Y = U ( L B );for i = 1 : length(Y) sprintf('%',Y(i)end end例8.5.1 取h 0.1 ,用差分法求下列边值问题的近似解y (x) 4(y(x) x), 0 x 1, y(0) 0, y(1) 2.解在MATLAB命令窗口输入：f = inline('-4.*x'); q= inline('4'); DiffMethod(f, q, 0, 1, , 0,

45、 2)回车，可的结果。8.5.2 用差分法解一般二阶常微分方程的边值问题%用差分法解一般二阶常微分方程的边值问题% a,b是区间,Step是步长,Alpha, Beta是初值% f, p, q 分别是式(8.2.12)中的 f (x), p (x), q (x)function Y = DiffMethod_2(f, p, q, a, b, Step, Alpha, Beta)N = (b-a) / Step;X = a : Step : b;A = zeros( N-1 );A(1, 1) = -2 * ( 2 - StepA2 * feval(p, X(2);A(1,2) = 2 + St

46、ep * feval(f, X(3);for i = 2 : N-2A(i,i) = -2 * ( 2 - feval(p, X(i+1) * StepA2 );A(i,i-1)= 2 - Step * feval(f, X(i+1);A(i-1,i) = 2 + Step * feval(f, X(i+1);endA(N-1,N-2) = 2 - Step * feval(f, X(N);A(N-1,N-1) = -2 * ( 2 - feval(p, X(N) * StepA2 );B(1) = 2*StepA2 * feval(q, X(2) - ( 2 - Step * feval(f

47、, X(2) ) * Alpha;B(2:N-2) = 2 * StepA2 * feval(q, X(3:N-1);B(N-1) = 2*StepA2 * feval(q, X(N) - ( 2 + Step * feval(f, X(N) ) * Beta;B = B'L,U = lu( A );Y = U ( L B );for i = 1 : length(Y) sprintf('%',Y(i)end end例8.5.2 取h /16 ,用差分法求下列边值问题的近似解y (x)y (x)2y(x) cosx, 0 x 2,y(0)0.3, y( 2)0.1.解在

48、MATLAB命令窗口输入：f = inline('cos(x)'); p = inline('-l'); q = inline('-2'); DiffMethod_2(f, p, q, 0, pi/ 2, pi/16,回车，可得结果。8.5.3 用有限元法解二阶常微分方程的边值问题%用有限元法解二阶常微分方程的边值问题%a,b是区间，h是步长，Alpha, Beta是初值 % f, p, q 分别是式(8.3.1)中的 f (x), p (x), q (x) function Y = FiniElem(f, p, q, a , b, Step,

49、Alpha, Beta) N = length(Step) - 1;X(1) = a;for i = 1 : N+1X(i+1) = X(i) + Step(i);endsyms t;ff = -1/Step(N+1) * feval(f, X(N)+t*Step(N+1) + . Step(N+1)* feval(p, X(N)+t*Step(N+1) * (1-t)*t;Knnn = int(ff, t, 0, 1); I syms t;ff = -1/Step(2) * feval(f, X(1)+t*Step(2) + . Step(2)*feval(p, X(1)+t*Step(2)*(1-t)*t;K101 = int(ff, t, 0, 1);for i = 2 : Nsyms t;ff = Step(i+1) * feval(q, X(i)+t*Step(i+1) * (1-t);B1(i) = int(ff, t, 0,1);endfor i = 1 : N-1syms t;ff = Step(i+1) * feval(q, X(i)+t*Step(i+1) * t;B2(i) = int(ff, t, 0, 1