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

摘    要 本文主要研究二阶常微分方程边值问的数值解法。对线性边值问题，我们总结了两类常用的数值，即打靶法和有限差分方法，对每种方法都列出了详细的计算步骤和Matlab程序代码，通过具体的算例对这两类方法的优缺点进行了细致的比较。 关键字：常微分方程边值问题；打靶法；差分法； ABSTRACT This article mainly discusses the numerical methods for solving Second-Order boundary value problems for Ordinary Differential Equations. On the one hand, we review two types of commonly used numerical methods for linear boundary value problems, i.e. shooting method and finite difference method. For each method, we give both the exact calculating steps , we compare the advantages and disadvantages in detail of these two methods through a specific numerical example. Key words：Boundary-Value Problems for Ordinary Differential Equations；Shooting Method；Finite Difference Method； 目    录 第一章  引言    - 1 - 第二章  二阶线性常微分方程    - 3 - 2.1试射法（“打靶”法）    - 4 - 2.1.1简单的试射法    - 4 - 2.1.2 基于叠加原理的试射法    - 5 - 2.2 有限差分法    - 11 - 2.2.1 有限差分逼近的相关概念    - 12 - 2.2.2 有限差分方程的建立    - 14 - 2.2.3 其他边值条件的有限差分方程    - 15 - 2.2.4 有限差分方程的解法    - 17 - 第三章  二阶非线性微分方程    - 22 - 3.1基于牛顿迭代法的打靶法    - 22 - 3.1.1 第一类边值条件推导    - 22 - 3.1.2  其他边值条件的推导    - 24 - 3.1.3  算法及程序代码    - 25 - 3.2 基于改进的牛顿迭代法的打靶法    - 31 - 3.2.1  算法的推导    - 31 - 3.2.2  算法及代码    - 32 - 第四章  改进算法的算例    - 38 - 第五章 总结    - 45 - 参考文献    - 46 - 致谢    - 47 - 第一章  引言 微分方程是现代数学中一个很重要的分支，从早期的微积分时代起，这个学科就成为了理论研究和实践应用的一个重要领域。在微分方程理论中，定解条件通常有两种提法：一种是给出了积分曲线在初始时刻的性态，相应的定解条件称为初值问题；另一种是给出了积分曲线首末两端的性态，这类条件则称为边界条件，相应的定解问题称为边值问题。 常微分方程边值问题在应用科学与工程技术中有着非常重要的应用，例如工程学、力学、天文学、经济学以及生物学等领域中的许多实际问题通常会归结为常微分方程边值问题[12]的求解。文献[9]给出了边值问题求解的方法，虽然求解常微分方程边值问题有很多解析方法可以求解，但这些方法只能用来求解一些特殊类型的方程，对从实际问题中提炼出来的微分方程往往不再适用，因而对常微分方程边值问题的数值方法的研究显得尤为重要。经典的数值方法主要有：试射法（打靶法）和有限差分法，见文献[2]。对于二阶线性边值问题，差分法的优点在于稳定性较好，但它的精度不高。而用打靶法求解线性问题时，解的精度较高，这是因为打靶法将边值问题的求解转化为相应的初值问题的求解，因而可以使用具有较高精度的Runge-Kutta法（见文献[1]），但是算法稳定性较差。 在本文中，我们首先总结了二阶线性边值问题的数值算法：打靶法、有限差分法。对每种方法都列出了详细的计算步骤和Matlab程序代码，通过具体的算例对这两类方法的优缺点进行了细致的比较。由于简单的打靶法过分依赖经验，我们考虑了基于线性叠加原理的打靶法，将线性边值问题转化为两个初值问题，并通过线性叠加得到原边值问题的解。 第二章  二阶线性常微分方程 二阶常微分方程一般可表示成如下的形式： ，                      (1) 边值条件有如下三类[9]： 第一类边值条件 ，                          (2) 第二类边值条件 ，                          (3) 第三类边值条件[19] ，                (4) 其中，  ，  ， 。 在对边值问题用数值方法求解之前，应该从理论上分析该边值问题的解是否存在，若问题的解不存在，用数值方法计算出来的数据没有任何意义。下面的定理给出了边值问题存在唯一解的充分条件。 定理1.1设方程(2.1)中的函数及，在区域 内连续，并且 (ⅰ)  ； (ⅱ) 在内有界，即存在常数，使得 ，  ， 则边值问题(2.1)-(2.4)的解存在且唯一[18]。 本章我们假设函数可以简单地表示成 ， 即边值问题(2.1)-(2.2)为具有如下形式的二阶线性边值问题             (2.5) 此时，解的存在唯一性定理(定理1.1)可以简单地表述为下述推论。 推论2.1  设，，在[a，b]上连续，且在[a，b]内，则线性边值问题(2.5)的解存在且唯一. 注：如无特别说明，本文中用到的(2.5)中的微分方程均满足推论2.1中的条件。     对于一般的边值问题(2.5)，尽管我们可以从理论上保证解的存在唯一性，但通常很难用解析方法求出其精确解，从实际问题中归结出来的微分方程往往不是解析方法所适用的特殊类型，此时数值方法显得尤为重要。下面我们先给出求解边值问题(2.5)的两种常用的数值方法，即试射法和有限差分法。 2.1试射法（“打靶”法） 2.1.1简单的试射法    对于边值问题(2.5)的求解，“试射法”的基本思想是将边值问题转化成初值问题来求解，即根据边界条件(2.2)，寻求与之等价的初始条件(2.6)，也就是说，反复是调整初始时刻的斜率值，使得初值问题(2.7)的积分曲线能“命中” . 设凭借经验能够提供的两个预测值，，我们按这两个斜率“试射”， 通过Runge-Kutta算法求解相应的初值问题(2.7)可以得到的两个预测值分别为，。若和都不满足预定的精度，则可用线性插值的方法校正，得到新的斜率值                       (2.8) 然后再按斜率试射，求解相应的初值问题(2.7)又得到新的结果.若或[16]，则可将作为的近似值；否则，继续过程(2.8)直到找到计算结果与相当符合为止。 上述试射法中初值的选取过分依赖于经验，局限性很大。下面介绍基于线性叠加原理的试射法。 2.1.2 基于叠加原理的试射法 设二阶线性常微分方程边值问题(2.5)的解存在并且唯一，并定义线性算子： .                      (2.9) 由于线性微分方程的解具有叠加性，即非齐次线性微分方程的解可由它的一个特解和相应的齐次线性微分方程的解的线性组合来表示，因此我们可以考虑如下的两个线性微分方程的初值问题：                           (2.10) 和                             (2.11) 事实上，设和分别为问题(2.10)和(2.11)的解，则不难验证                       (2.12) 是问题(2.5)的解，其中. 通过上述描述，我们可以得到基于叠加原理的打靶法的基本步骤为： 1. 根据边值问题(2.5)构造相应的初值问题(2.10)和(2.11)； 2. 分别求出两个初值问题(2.10)和(2.11)的解和； 3. 将和按(2.12)式做组合，所得的函数就是边值问题(2.5)的解.其几何图像见图2.1。 图 2.1  打靶法求解结果的几何图像 (2.10)和(2.11)均为二阶常微分方程初值问题，求解时可通过引入变量代换将其化成相应的一阶方程组初值问题。如令：                             (2.13) 则(2.10)式可以写成                       (2.14) (2.11)式可以写成                           (2.15) 这样就可以利用Runge-Kutta方法求解(2.14)和(2.15)。 对于更一般的线性边值问题：             (2.16) 用基于叠加原理的打靶法的步骤为： 1. 根据(2.16)式，构造两个相应的初值问题：                           (2.17) 和                               (2.18) 其中和是满足条件的两个任意的常量. 2. 求解初值问题(2.17)和(2.18)式，设其解分别为和. 3. 将和做线性组合 由此计算得到的函数就是(2.16)式的解，其具体的数值计算格式可描述为如下算法： 算法1  基于叠加原理的试射法求解线性边值问题(2.5)的算法如下： 1. 输入 区间的端点，，及边界条件常数，，区间等分数 2. 计算 步长，，，， 3.  for i=0，1，2，…，N-1     3.1     3.2 计算             3.3 计算        end for(i) 4. 计算 5. 输出，， 6. for i=1，2，…，N     6.1 计算     6.2 输出，，(说明：，分别是，的近似值)   end for(i) 代码1 将上面的算法在matlab实现代码如下： function [x,y]=lsh(func1,func2,a,b,ya,yb,N) %利用RK方法计算二阶线性微分方程的边值问题， %由微分方程知识可以知道，对于二阶微分方程可以化为一阶微分方程组 %func1为第一个微分方程组第二个函数，func2为第二个微分方程组第二个函数 %初始时刻为a，结束时为b，初始端点函数值为ya，末端点函数值yb %N为分成的区间数目 h=(b-a)/N; u(1,1)=ya; u(1,2)=0; v(1,1)=0; v(1,2)=1; for i=1:N     x(i,:)=a+(i-1)*h;     K1=h*u(i,2);     L1=h*feval(func1,x(i,:),u(i,1),u(i,2));     K2=h*(u(i,2)+1/2*L1);     L2=h*feval(func1,x(i,:)+1/2*h,u(i,1)+1/2*K1,u(i,2)+1/2*L1);     K3=h*(u(i,2)+1/2*L2);     L3=h*feval(func1,x(i,:)+1/2*h,u(i,1)+1/2*K2,u(i,2)+1/2*L2);     K4=h*(u(i,2)+L3);     L4=h*feval(func1,x(i,:)+h,u(i,1)+K3,u(i,2)+L3);     u(i+1,1)=u(i,1)+1/6*(K1+2*K2+2*K3+K4);     u(i+1,2)=u(i,2)+1/6*(L1+2*L2+2*L3+L4);         k1=h*v(i,2);     l1=h*feval(func2,x(i,:),v(i,1),v(i,2));     k2=h*(v(i,2)+1/2*l1);     l2=h*feval(func2,x(i,:)+1/2*h,v(i,1)+1/2*k1,v(i,2)+1/2*l1);     k3=h*(v(i,2)+1/2*l2);     l3=h*feval(func2,x(i,:)+1/2*h,v(i,1)+1/2*k2,v(i,2)+1/2*l2);     k4=h*(v(i,2)+l3);     l4=h*feval(func2,x(i,:)+h,v(i,1)+k3,v(i,2)+l3);     v(i+1,1)=v(i,1)+1/6*(k1+2*k2+2*k3+k4);     v(i+1,2)=v(i,2)+1/6*(l1+2*l2+2*l3+l4);    end u v x(N+1,:)=x(N,:)+h; y(1,1)=ya; y(1,2)=(yb-u(N+1,1))/v(N+1,1); for i=1:N     y(i+1,1)=u(i+1,1)+y(1,2)*v(i+1,1);     y(i+1,2)=u(i+1,2)+y(1,2)*v(i+1,2); end 例2.1 用基于叠加原理的打靶法求解下面的线性边值问题： 解：该边值问题有精确解[8,21] 其中，. 对该问题用算法1求解时需要近似求解初值问题 和 当，时，计算结果如下表所示。表中列出的值逼近的值，逼近的值，且值逼近 ， 表示逐点绝对误差。 表 2.1  基于叠加原理的打靶法算例数据 () 1.0 1.000000000 0.000000000 1.000000000 1.000000000 0.000000000 1.1 1.008960575 0.091179858 1.092629164 1.092629298 -0.134347735 1.2 1.032454720 0.168511750 1.187084707 1.187084840 -0.133672730 1.3 1.066743750 0.236087037 1.283382266 1.283382364 -0.097795783 1.4 1.109287947 0.296590670 1.381445892 1.381445951 -0.060163483 1.5 1.158299996 0.351843791 1.481159386 1.481159416 -0.030633643 1.6 1.212483715 0.403116946 1.582392450 1.582392460 -0.010770011 1.7 1.270874543 0.451318399 1.685013962 1.685013961 0.000543515 1.8 1.332738514 0.497111366 1.788898540 1.788898534 0.005050135 1.9 1.397506177 0.540989278 1.893929514 1.893929509 0.004410489 2.0 1.464728147 0.583325383 2.000000000 2.000000000 0.000000000 其精确解和数值解的比较如下图所示 图2.2  基于叠加原理的打靶法算例图像 图像中的直线表示边值问题的精确解，蓝色点表示利用算法1得到的数值解 2.2 有限差分法 有限差分方法是用于微分方程定解问题求解的最广泛的数值方法，其基本思想是用离散的、只含有有限个未知量的差分方程去近似代替连续变量的微分方程和定解条件，并把相应的差分方程的解作为微分方程定解问题的近似解。 本节我们依然讨论边值问题(2.5)，介绍解两点边值问题的另一种数值方法——有限差分方法。 2.2.1 有限差分逼近的相关概念 设函数光滑，且，利用Taylor展开，可得             (2.19)             (2.20) 由(2.19)可以得到一阶导数的表达式           (2.21a) 或者                     (2.21b) 同理由(2.20)式可得             (2.22a) 或者                     (2.22b) 其中表示截断误差项.因此，可得一阶导数的的差分近似表达式为                         (2.23)                         (2.24) 由(2.21)和(2.22)可知，差商(2.23)和(2.24)逼近微商的精度为一阶，即为，为了得到更精确的差分表达式，将(2.19)减(2.20)可得               (2.25) 从而可以的到                 (2.26a) 或者                   (2.26b) 其中，. 可得一阶导数的差分近似表达式为                       (2.27) 由此可知，(2.16)差商逼近微商的精度为二阶，即为。 公式(2.23)，(2.24)和(2.27)分别被称为逼近一阶微商的向前，向后和中心差分公式[6]。 类似地，我们还可以给出二阶微商和高阶微商的差分近似表达式。例如将(2.19)和(2.20)两式相加可得 进而有               (2.28) 其中. 因此，二阶导数的差分近似表达式[8]为                   (2.29) 下表中列出了一些高阶导数的常用差分近似公式以及其相应的截断误差阶： 表 2.2  导数的差分逼近及其相应的截断误差阶 导数(或微商) 有限差分逼近[9,22] 截断误 差阶           2.2.2 有限差分方程的建立 考虑边值问题(2.5)，取一正整数将区间[a，b]分成等份。设其节点分别为，。其中为步长，为内节点，，称为边界节点。在每一个节点处，边值问题的方程(2.5)可以转化为                     (2.30) 假设，利用上一节内容可得到和的中心差商公式：                     (2.31) 其中，.               (2.32a) 或者           (2.32b) 其中，. 将(2.27)和(2.32a)代入(2.30)式可得 (2.33) 其中，                     (2.34) 称为用差分方程逼近微分方程的截断误差[9]。略去，并用代替可得到逼近边值问题(2.30)式的差分方程：             (2.35) 再记                       (2.36) 则由(2.35)整理可得到关于的方程组：       (2.37) 该方程组的可解性和解的误差估计可参考文献[10,20]。 2.2.3 其他边值条件的有限差分方程 在实际的应用中，我们可能会碰到边值条件为(2.3)或(2.4)的边值问题，即第二或第三边值条件.对于第二或第三边值条件，在建立差分方程时，必须选用适当的差商公式代替边界上的和。例如，为了使近似解在边界和上具有二阶精度，可以选用本章第一节表2.2中列出的差商公式近似：             (2.38) 在(2.31)中，取近似 那么我们可以得到第二边值条件的差分方程组:         (2.39) 第三边值条件的差分方程组：         (2.40) 其中包含个未知量。其系数矩阵都具有下列形态. 其中，*代表不恒为零的元素。由线性代数的知识可以知道，矩阵可经过适当的初等行变换可化为三对角矩阵。 2.2.4 有限差分方程的解法 下面我们来研究下这些差分方程组的解法（以第一边值条件所得的差分方程组(2.37)为例）： 将(2.37)简记为                               (2.41) 其中， ， 从中可以看出为三对角矩阵，可以用追赶法[14]求解。 算法2  用有限差分法求解二阶线性边值问题(2.5)的算法如下： 1. 输入区间的端点，，及边界条件常数，，区间等分数 2. 计算（形成方程组(2.41)中的矩阵，及右端的向量，为节省空间，用代替）                 for i=2:N-2                     end for(i)         3.  用追赶法计算三对角方程组(2.41)中的矩阵[2,3] 3.1 计算 for i=2，3，…，N-2         end for(i)         3.2 解 （该存储的值）     for  i=2，3，…，N-1         end for(i) 3.3 解（为了节省空间用继续用存储(2.41)的解）     for i=N-2，…，2，1         end for(i) 4  输出相应的值. 例2.2  用有限差分法解例题2.1中的边值问题： 解：考虑有限差分法来求解这个问题。在算法2中，取，。计算结果见下表2.3及其图像见图2.3，值逼近，表示逐点绝对误差。 表 2.3  有限差分法算例数据 () 1.0 1.000000000 1.000000000 0.000000000 1.1 1.092600520 1.092629298 -0.287777612 1.2 1.187043128 1.187084840 -0.417116886 1.3 1.283336870 1.283382364 -0.454938647 1.4 1.381402046 1.381445951 -0.439054684 1.5 1.481120262 1.481159417 -0.391548876 1.6 1.582359895 1.582392461 -0.325650629 1.7 1.684989018 1.685013962 -0.249433495 1.8 1.788881746 1.788898535 -0.167884482 1.9 1.893921099 1.893929509 -0.084100540 2.0 2.000000000 2.000000000 0.000000000 精确解与用有限差分法求得的数值解的比较如下图所示 图 2.3  有限差分法算例图像 图像中的直线表示边值问题的精确解，蓝色点表示由算法2得到的数值解. 从该题的结果中我们可以看出，例2.2的数值解的精度与例2.1相比差很多。原因是在例2.1中所用的方法是建立在Runge-Kutta算法的基础上的，而Runge-Kutta算法具有的精度，而有限差分法的局部截断误差阶仅有。 第五章 总结 通过本文的研究，我们总结了求解二阶线性常微分方程边值问题以及二阶线性边值问题的数值解法。对于二阶线性边值问题时，本文总结了三种方法，并且得出结论：简单的“打靶法”过分依赖于经验，盲目性很大；基于“叠加原理”的打靶法是建立在经典的Runge-Kutta算法的基础上，因此在求解线性边值问题时具有较好的精度，但是稳定性相对较差；有限差分法，可以适应很多线性边值问题的求解，因此稳定性较好，但是其精度较基于“叠加原理”的打靶法要差些。 参 考 文 献 [1] 李庆扬等. 数值分析 (第四版)[M].北京：清华大学出版社，2001.08； [2] Richard L.Burden  J.Douglas Faires.  NUMERICAL ANALYSIS (Seventh Edition)[M]. Beijing：Higher Education, 2001.08； [3] Rainer Kress. Numerical Analysis [M]. New York: Springer-verlag, 1998； [4] David Kincaid. Numerical Analysis: Mathematics of Scientific Computing[M].Beijing：China Machine Press，2003.4; [5] H.J.Pech. Multipoint Boundary-Value Solutions of Two-Point Boundary-Value Problems[M]. JOURNAL OF OPTIMIZATION THEORY AND APPLICATIONS. 1999.01； [6] 李庆扬等. 数值分析 (第三版)[M].北京：清华大学出版社，1986.12； [7] 黄争鸣. 两点边值问题的一种有效算法[J]. 计算物理，1992.06； [8] 李磊. 边值问题的并行算法[J]. 工程数学学报，1989.03； [9] 葛渭高等. 常微分方程与边值问题[M]. 北京：科学出版社，2008； [10] 余德浩. 微分方程数值解法[M]. 北京：科学出版社，2003.10； [11] John H.Mathews Kurtis D.Fink Numerical Methods Using MATLAB Fourth Edition[M]. Beijing：Publishing House of Electronics Industry. [12] 王高雄等.常微分方程第三版[M].北京：清华大学出版社，2006； [13] 宋叶志等.MATLAB数值分析与应用[M].北京：机械工业出版社，2009； [14] Curtis F.Gerald Applied Numerical Analysis(Seventh Edition)[M]. Beijing：Higher Education Press,2005.05； [15] 刘墨德.  Newton 迭代法及其改进[J]. 三明学院学报，2007.06； [16] 赵秋玲等. 求解二点边值问题打靶法的一种改进方法[J]. 上海力学，1999.12 [17] 李丽容. 对牛顿迭代法的改进[J]. 中国水运(理论版)，2006.05； [18] 张春梅. 一类二阶非线性常微分方程边值问题的有效解法[J]. 优秀硕士，2006.06； [19] 徐静等. 非线性一般两点边值问题的试射法[J]. 内江科技，2006.08； [20] 帕力旦·赛力提尼亚孜等，试射法在求解二阶线性微分方程中的应用[J]. 新疆大学学报(自然科学版)，2006.02； [21] 向隆万. 两点边值问题的双向迭加法[J]. 工程数学学报，1984.04 [22] 周保民. 常微分方程边值问题的高精度差分法[J].