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

第8章 常微分方程边值问题的数值解法8.1 引 言第7章介绍了求解常微分方程初值问题的常用的数值方法；本章将介绍常微分方程的边值问题的数值方法。只含边界条件(boundary-value condition)作为定解条件的常微分方程求解问题称为常微分方程的边值问题(boundary-value problem). 为简明起见，我们以二阶边值问题为例介绍常用的数值方法。一般的二阶常微分方程边值问题(boundary-value problems for second-order ordinary differential equations)为, (8.1.1)其边界条件为下列三种情况之一：(1) 第一类边界条件 (the first-type boundary conditions)：(2) 第二类边界条件 (the second-type boundary conditions)：(3) 第三类边界条件 (the third-type boundary conditions)： 定理8.1.1 设(8.1.1)中的函数及其偏导数, 在上连续. 若(1) 对所有，有；(2) 存在常数，对所有，有,则边值问题(8.1.1)有唯一解。推论 若线性边值问题 (8.1.2)满足(1) 和在上连续；(2) 在上, ,则边值问题(8.1.1)有唯一解。求边值问题的近似解，有三类基本方法：(1) 差分法(difference method)，也就是用差商代替微分方程及边界条件中的导数，最终化为代数方程求解；(2) 有限元法(finite element method)；(3) 把边值问题转化为初值问题，然后用求初值问题的方法求解。8.2 差分法8.2.1 一类特殊类型二阶线性常微分方程的边值问题的差分法设二阶线性常微分方程的边值问题为其中在上连续，且.用差分法解微分方程边值问题的过程是：(i) 把求解区间分成若干个等距或不等距的小区间，称之为单元；(ii) 构造逼近微分方程边值问题的差分格式. 构造差分格式的方法有差分法, 积分插值法及变分插值法；本节采用差分法构造差分格式；(iii) 讨论差分解存在的唯一性、收敛性及稳定性；最后求解差分方程.现在来建立相应于二阶线性常微分方程的边值问题(8.2.1), (8.2.2)的差分方程.( i ) 把区间等分，即得到区间的一个网格剖分：,其中分点，并称之为网格节点(grid nodes)；步长.( ii ) 将二阶常微分方程(8.2.2)在节点处离散化：在内部节点处用数值微分公式 (8.2.3)代替方程(8.2.2)中，得, (8.2.4)其中.当充分小时，略去式(8.2.4)中的，便得到方程(8.2.1)的近似方程, (8.2.5)其中, 分别是的近似值, 称式(8.2.5)为差分方程(difference equation)，而称为差分方程(8.2.5)逼近方程(8.2.2)的截断误差(truncation error). 边界条件(8.7.2)写成 (8.2.6)于是方程(8.2.5), (8.2.6)合在一起就是关于个未知量，以及个方程式的线性方程组： (8.2.7)这个方程组就称为逼近边值问题(8.2.1), (8.2.2)的差分方程组(system of difference equations)或差分格式(difference scheme)，写成矩阵形式. (8.2.8)用第2章介绍的解三对角方程组的追赶法求解差分方程组(8.2.7)或(8.2.8), 其解称为边值问题(8.2.1), (8.2.2)的差分解(difference solution). 由于(8.2.5)是用二阶中心差商代替方程(8.2.1)中的二阶微商得到的，所以也称式(8.2.7)为中心差分格式(centered-difference scheme).( iii ) 讨论差分方程组(8.2.7)或(8.2.8)的解是否收敛到边值问题(8.2.1), (8.2.2)的解，估计误差.对于差分方程组(8.2.7)，我们自然关心它是否有唯一解；此外，当网格无限加密，或当时，差分解是否收敛到微分方程的解. 为此介绍下列极值原理：定理8.2.1 (极值原理) 设是给定的一组不全相等的数，设 . (8.2.9)(1) 若, 则中非负的最大值只能是或；(2) 若, 则中非正的最小值只能是或.证 只证(1) 的情形，而(2) 的情形可类似证明. 用反证法. 记，假设, 且在中达到. 因为不全相等，所以总可以找到某个，使，而和中至少有一个是小于的. 此时因为，所以, 这与假设矛盾，故只能是或. 证毕！推论 差分方程组(8.2.7)或(8.2.8)的解存在且唯一.证明 只要证明齐次方程组 (8.2.10)只有零解就可以了. 由定理8.7.1知，上述齐次方程组的解的非负的最大值和非正的最小值只能是或. 而，于是 证毕！利用定理8.2.1还可以证明差分解的收敛性及误差估计. 这里只给出结果：定理8.2.2 设是差分方程组(8.2.7)的解，而是边值问题(8.2.1), (8.2.2)的解在上的值，其中. 则有 (8.2.11)其中.显然当时，. 这表明当时，差分方程组(8.2.7)或(8.2.8)的解收敛到原边值问题(8.7.1), (8.7.2)的解.例8.2.1 取步长，用差分法解边值问题 并将结果与精确解进行比较.解 因为，, 由式(8.2.7)得差分格式 , , 其结果列于表8.2.1.表8.2.1准确值010010.10. 03329230.033365620.20. 06491630.065060430.30. 09313690.093346140.40. 11608310.116348250.50. 13167250.131979660.60. 13752880.137857870.70. 13088630.131208780.80. 10847930.108755390.90. 06641140.0665865101.000从表8.2.1可以看出, 差分方法的计算结果的精度还是比较高的. 若要得到更精确的数值解，可用缩小步长的方法来实现.8.2.2 一般二阶线性常微分方程边值问题的差分法对一般的二阶微分方程边值问题 (8.2.12)假定其解存在唯一. 为求解的近似值，类似于前面的做法，( i ) 把区间等分，即得到区间的一个网格剖分：,其中分点，步长.( ii ) 对式(8.2.12)中的二阶导数仍用数值微分公式代替，而对一阶导数，为了保证略去的逼近误差为，则用3点数值微分公式；另外为了保证内插，在2个端点所用的3点数值微分公式与内网格点所用的公式不同，即 (8.2.13)略去误差，并用的近似值代替，便得到差分方程组 (8.2.14)其中, 是的近似值. 整理得 (8.2.15)解差分方程组(8.2.15)，便得边值问题(8.2.12)的差分解.特别地, 若，则式(8.2.12)中的边界条件是第一类边值条件：此时方程组(7.7.16)为 (8.2.16)方程组(8.2.16)是三对角方程组，用第2章介绍的解三对角方程组的追赶法求解差分方程组(8.2.16)，便得边值问题(8.2.12)的差分解.( iii ) 讨论差分方程组(8.2.16)的解是否收敛到微分方程的解，估计误差. 这里就不再详细介绍.例8.2.2 取步长，用差分法求下列边值问题的近似解，并将结果与精确解进行比较.精确解是.解 因为，, 由式(8.2.17)得差分格式 , , 其结果列于表8.2.2.表8.2.2准确值00-0.3-0.31/16-0.3137967-0.313744622/16-0.3154982-0.315432233/16-0.3050494-0.304997944/16-0.2828621-0.282842755/16-0.2497999-0.249818066/16-0.2071465-0.207193077/16-0.1565577-0.15660568/2-0.1000000-0.10000008.3 有限元法有限元法(finite element method)是求解微分方程定解问题的有效方法之一，它特别适用在几何、物理上比较复杂的问题. 有限元法首先成功地应用于结构力学和固体力学，以后又应用于流体力学、物理学和其他工程科学. 为简明起见，本节以线性两点边值问题为例介绍有限元法.考虑线性两点边值问题其中, . 此微分方程描述了长度为的可变交叉截面(表示为)的横梁在应力和下的偏差.8.3.1 等价性定理 记, 引进积分. (8.3.3)任取，就有一个积分值与之对应，因此是一个泛函(functional)，即函数的函数. 因为这里是的二次函数，因此称为二次泛函.对泛函(8.3.3)有如下变分问题(variation problem)：求函数，使得对任意, 均有, (8.3.4)即在处达到极小, 并称为变分问题(8.3.4)的解.可以证明：定理8.3.1(等价性定理) 是边值问题(8.3.1), (8.3.2)的解的充分必要条件是使泛函在上达到极小，即是变分问题(8.3.4)在上的解.证 (充分性) 设是变分问题的解；即使泛函在上达到极小，证明必是边值问题(8.3.1), (8.3.2)的解.设是任意一个满足的函数，则函数，其中为参数. 因为使得达到极小，所以，即积分作为的函数，在处取极小值，故. (8.3.5)计算上式，得利用分部积分法计算积分代入式(8.3.6)，得因为是任意函数，所以必有. (8.3.8)否则，若在上某点处有，不妨设，则由函数的连续性知，在包含的某一区间上有.作显然，且，但，这与式(8.3.7)矛盾. 于是式(8.3.8)成立，即变分问题(8.3.4)的解满足微分方程(8.3.1), 且故它是边值问题(8.3.1), (8.3.2)的解.(必要性) 设是边值问题(8.3.1), (8.3.2)的解，证明是变分问题(8.3.4)的解；即：使泛函在上达到极小.因为满足方程(8.3.1)，所以. 设任意，则函数满足条件，且. 于是因为，所以当时，, 即.只有当时，. 这说明使泛函在上达到极小. 证毕！定理8.3.2 边值问题(8.3.1), (8.3.2)存在唯一解.证明 用反证法. 若都是边值问题(8.3.1), (8.3.2)的解，且不相等，则由定理8.3.1知，它们都使泛函在上达到极小，因而 且 ，矛盾！因此边值问题(8.3.1), (8.3.2)的解是唯一的. 由边值问题解的唯一性，不难推出边值问题(8.3.1), (8.3.2)解的存在性(留给读者自行推导). 8.3.2 有限元法等价性定理说明，边值问题(8.3.1), (8.3.2)的解可化为变分问题(8.3.4)的求解问题. 有限元法就是求变分问题近似解的一种有效方法. 下面给出其解题过程：第1步 对求解区间进行网格剖分区间称为单元，长度称为步长，. 若，则称上述网格剖分为均匀剖分. 给定剖分后，泛函(8.3.3)可以写成. (8.3.9)第2步 构造试探函数空间。为了计算积分(8.3.9)，最简单的近似方法是将分段线性函数的集合作为试探函数空间。设分别是边值问题(8.3.1), (8.3.2)的解在节点处的值，用分段线性插值 (8.3.10)近似，并称式(8.3.10)为单元形状函数(element shape function).为了将线性插值(8.3.10)标准化，令，则将变到轴上的单元. 记，则式(8.3.10)可写成 (8.3.11)第3步 建立有限元方程组. 将式(8.3.10)代入泛函(8.3.9)，有.由式(8.3.11)知 (8.3.12)式(8.3.12)右端第1个求和号内的第项(对应第个单元)是关于的二次型，可以写成, (8.3.13)其中，, (8.3.14)称为单元刚度矩阵(element stiffness matrix)， (8.3.15)而式(8.3.12)的第2个求和号内的第项可以写成 (8.3.16)其中，于是式(8.3.12)中求和号内的项可以写成. (8.3.17)再令及矩阵则. 于是式(8.3.17)又可写成 . (8.3.18)把式(8.3.18)代入式(8.3.12)右端求和号内，得. (8.3.19)记，(8.3.20), (8.3.21)则式(8.3.19)化为 (8.3.22)并称为总刚度矩阵(total stiffness matrix)，称为右端向量. 因为是使取极小值的函数，所求的自然使式(8.3.22)的右边取极小值，所以应有. (8.3.23)记, 则式(8.8.23)为或 (8.3.24)得方程组. (8.3.25)因为是已知的，不能任意选取，所以不能要求式(8.3.23)对也成立. 因此方程组(8.3.24)或(8.3.25)中应当去掉首末2个方程，并把其他方程中含有的改为已知量，所得方程组就是未知量满足的代数方程组. (8.3.26)方程组(8.3.25)或(8.3.26)称为有限元方程组(finite element system of equations). 用第2章介绍的解三对角方程组的追赶法求解有限元方程组(8.3.26)，可解出，即得变分问题(8.3.4)的近似解，也就是边值问题(8.3.1), (8.3.2)解的近似值. 这种方法称为有限单元法(finite element method)或有限元法.例8.3.1 用有限元法求下列边值问题的近似解，并将结果与精确解进行比较.取，精确解是.解 因为，, 由式(8.3.26)得有限元方程组 其结果列于表8.3.1.表8.3.1准确值000010.2-0.1644-0.1620.4-0.2443-0.2430.6-0.2434-0.2440.8-0.1620-0.165100上面虽然是对边值问题(8.3.1), (8.3.2)介绍的有限元解法，但其解题步骤对一般的微分方程定解问题也是适用的. 对所给微分方程定解问题，首先找出相应微分方程的变分形式，然后进行下列步骤：第1步 对定义区域(或定义区间)进行网格剖分, 将定义区域(或定义区间)剖分为若干个小单元(一维是区间；多维是区域，如矩形、三角形等)；第2步 构造试探空间; 即选择适当的插值函数类.第3步 建立有限元方程组. 计算单元刚度矩阵及右端向量，再形成总刚度矩阵及总右端项，写出有限元方程组；结合定解条件，求解有限元方程组.注 从形式上看，用有限元法解微分方程定解问题很繁，但是有限元法的求解步骤规范，便于在计算机上实现，并且总刚度矩阵是三对角对称正定矩阵，因此有限元方程组可用第2章介绍的解三对角方程组的追赶法求解. 有限元法最主要的优点是可以处理相当复杂的区域上的边值问题。8.4 打靶法解常微分方程边值问题的打靶法(shooting method)，也称为尝试法，其基本思想是把边值问题转化为初值问题来解：首先作出一些只满足一端边值条件的解，然后再从这些解中找出适合另一端边值条件的解. 下面以二阶常微分方程带第一类边界条件的边值问题为例介绍常微分方程边值问题的打靶法.7.6节曾介绍过二阶常微分方程初值问题(7.6.10) 的求解方法. 将上式中的改为，改为，得 (8.4.3)设初值问题(8.4.3)的解为，显然依赖于，即. 为了求解边值问题(8.4.1), (8.4.2)，必须求，使之满足. 下面介绍2种方法来求.方法1 根据实际问题情况选一个，解初值问题 (8.4.4)得到的解仍记为. 若或(为事先给定的精度)，则就是边值问题(8.4.1), (8.4.2)的解. 否则，根据与的误差，将修改为；例如取, 再解初值问题 (8.4.5)得到解. 若满足或，则就是边值问题(8.4.1), (8.4.2)的解. 否则，再将适当修改为; 例如可用作线性插值求：,然后解初值问题 (8.4.6)的解. 如此继续下去，直到达到精度要求为止.方法2 求另一种方法是求函数的一个零点. 对于每一个自变量，通过解初值问题(8.4.4)，可解出. 计算, (8.4.7)然后采用第3章介绍的求方程根的方法求的零点. 首先寻找，使，则在区间或内至少有的一个零点. 然后可用二分法求. 也可用Newton迭代公式求的近似值.几何解释：微分方程边值问题(8.4.1), (8.4.2)的解是一条通过两点的曲线(如图8.4.1). 初值问题(8.4.4)的解是一条通过点、斜率为的曲线(如图8.4.1). 初值问题(8.4.5)的解是一条通过点、斜率为的曲线(如图8.4.1). 这有点象射击者从定点向目标射击一样. 根据经验以某一角度(斜率为)试射一次. 如果与目标相差太大，调整射击角度(斜率为)，再进行射击. 如此继续进行下去，直到击中目标，或击中的点与的误差在允许的范围之内。图8.4.1参考文献1提供的资料还讨论了选择初始值的重要性，并介绍了多重打靶法，这里就不作详细介绍，有兴趣的读者可参看相关资料.8.5 算法程序8.5.1 用差分法解二阶线性常微分方程的边值问题%用差分法解二阶线性常微分方程的边值问题%a,b是区间, Step是步长, Alpha, Beta是初值% f, q分别是式(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-1 A(i,i) = -1 * ( 2+feval(q, X(i+1) * Step2 ); if i = N-1 A(i,i+1) = 1; A(i+1,i) = 1; endendB(1) = Step2 * feval(f, X(2) - Alpha;B(2:N-2) = Step2 * feval(f, X(3:N-1);B(N-1) = Step2 * feval(f, X(N) - Beta;B = B;L,U = lu( A );Y = U ( L B ) ;for i = 1 : length(Y) sprintf(%10.7f,Y(i)endend例8.5.1 取，用差分法求下列边值问题的近似解解 在MATLAB命令窗口输入：f = inline(-4.*x); q= inline(4); DiffMethod(f, q, 0, 1, 0.1, 0, 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 - Step2 * feval(p, X(2) );A(1,2) = 2 + Step * feval(f, X(3);for i = 2 : N-2 A(i,i) = -2 * ( 2 - feval(p, X(i+1) * Step2 ); 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) * Step2 );B(1) = 2*Step2 * feval(q, X(2) - ( 2 - Step * feval(f, X(2) ) * Alpha;B(2:N-2) = 2 * Step2 * feval(q, X(3:N-1);B(N-1) = 2*Step2 * 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(%10.7f,Y(i)endend例8.5.2 取，用差分法求下列边值问题的近似解解 在MATLAB命令窗口输入：f = inline(cos(x); p = inline(-1); q = inline(-2);DiffMethod_2(f, p, q, 0, pi/2, pi/16, -0.3, -0.1)回车，可得结果。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, Alpha, Beta)N = length(Step) - 1;X(1) = a;for i = 1 : N+1 X(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);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 : N syms 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-1 syms t; ff = Step(i+1) * feval(q, X(i)+t*Step(i+1) * t; B2(i) = int(ff, t, 0, 1);endB(1) = B2(1) + B1(2) - Alpha*K101;B(N-1) = B2(N-1) + B1(N) - Beta*Knnn;for i = 2: N-2 B(i) = B2(i) + B1(i+1);e