1、微分方程是现代数学中一个很重要的分支，从早期的微积分时代起，这个学科就成为了理论研究和实践应用的一个重要领域。在微分方程理论中，定解条件通常有两种提法：一种是给出了积分曲线在初始时刻的性态，相应的定解条件称为初值问题；另一种是给出了积分曲线首末两端的性态，这类条件则称为边界条件，相应的定解问题称为边值问题。常微分方程边值问题在应用科学与工程技术中有着非常重要的应用，例如工程学、力学、天文学、经济学以及生物学等领域中的许多实际问题通常会归结为常微分方程边值问题12的求解。文献9给出了边值问题求解的方法，虽然求解常微分方程边值问题有很多解析方法可以求解，但这些方法只能用来求解一些特殊类型的方程，对
2、从实际问题中提炼出来的微分方程往往不再适用，因而对常微分方程边值问题的数值方法的研究显得尤为重要。经典的数值方法主要有：试射法（打靶法）和有限差分法，见文献2。对于二阶线性边值问题，差分法的优点在于稳定性较好，但它的精度不高。而用打靶法求解线性问题时，解的精度较高，这是因为打靶法将边值问题的求解转化为相应的初值问题的求解，因而可以使用具有较高精度的Runge-Kutta法（见文献1），但是算法稳定性较差。在本文中，我们首先总结了二阶线性边值问题的数值算法：打靶法、有限差分法。对每种方法都列出了详细的计算步骤和Matlab程序代码，通过具体的算例对这两类方法的优缺点进行了细致的比较。由于简单的打
3、靶法过分依赖经验，我们考虑了基于线性叠加原理的打靶法，将线性边值问题转化为两个初值问题，并通过线性叠加得到原边值问题的解。第二章 二阶线性常微分方程二阶常微分方程一般可表示成如下的形式：， （1）边值条件有如下三类9：第一类边值条件， （2）第二类边值条件， （3）第三类边值条件19， （4）其中， ， ， 。在对边值问题用数值方法求解之前，应该从理论上分析该边值问题的解是否存在，若问题的解不存在，用数值方法计算出来的数据没有任何意义。下面的定理给出了边值问题存在唯一解的充分条件。定理1.1设方程（2.1）中的函数及，在区域内连续，并且（） ；（） 在内有界，即存在常数，使得， ，则边值问题（
4、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），尽管我们可以从理论上保证解的存在唯一性，但通常很难用解析方法求出其精确解，从实际问题中归结出来的微分方程往往不是解析方法所适用的特殊类型，此时数值方法显得尤为重要。下面我们先给出求解边值问题
5、（2.5）的两种常用的数值方法，即试射法和有限差分法。2.1试射法（“打靶”法）2.1.1简单的试射法 对于边值问题（2.5）的求解，“试射法”的基本思想是将边值问题转化成初值问题来求解，即根据边界条件（2.2），寻求与之等价的初始条件（2.6），也就是说，反复是调整初始时刻的斜率值，使得初值问题（2.7）的积分曲线能“命中” .设凭借经验能够提供的两个预测值，我们按这两个斜率“试射”， 通过Runge-Kutta算法求解相应的初值问题（2.7）可以得到的两个预测值分别为，。若和都不满足预定的精度，则可用线性插值的方法校正，得到新的斜率值 （2.8）然后再按斜率试射，求解相应的初值问题（2.7
6、）又得到新的结果.若或16，则可将作为的近似值；否则，继续过程（2.8）直到找到计算结果与相当符合为止。上述试射法中初值的选取过分依赖于经验，局限性很大。下面介绍基于线性叠加原理的试射法。2.1.2 基于叠加原理的试射法设二阶线性常微分方程边值问题（2.5）的解存在并且唯一，并定义线性算子：. （2.9）由于线性微分方程的解具有叠加性，即非齐次线性微分方程的解可由它的一个特解和相应的齐次线性微分方程的解的线性组合来表示，因此我们可以考虑如下的两个线性微分方程的初值问题： （2.10）和 （2.11）事实上，设和分别为问题（2.10）和（2.11）的解，则不难验证 （2.12）是问题（2.5）的
7、解，其中.通过上述描述，我们可以得到基于叠加原理的打靶法的基本步骤为：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）和
8、（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. 计算
9、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
10、;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（
11、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（
12、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）;enduvx（N+1,:）=x（N,:）+h;y（1,1）=ya;y（1,2）=（yb-u（N+1,1）/v（N+1,1）; 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）;例2.1 用基于叠加原
13、理的打靶法求解下面的线性边值问题：解：该边值问题有精确解8,21其中，.对该问题用算法1求解时需要近似求解初值问题当，时，计算结果如下表所示。表中列出的值逼近的值，逼近的值，且值逼近表示逐点绝对误差。表 2.1 基于叠加原理的打靶法算例数据（）1.01.0000000000.0000000001.11.0089605750.0911798581.0926291641.092629298-0.1343477351.21.0324547200.1685117501.1870847071.187084840-0.1336727301.31.0667437500.2360870371.28338226
14、61.283382364-0.0977957831.41.1092879470.2965906701.3814458921.381445951-0.0601634831.51.1582999960.3518437911.4811593861.481159416-0.0306336431.61.2124837150.4031169461.5823924501.582392460-0.0107700111.71.2708745430.4513183991.6850139621.6850139610.0005435151.81.3327385140.4971113661.7888985401.788
15、8985340.0050501351.91.3975061770.5409892781.8939295141.8939295090.0044104892.01.4647281470.5833253832.000000000其精确解和数值解的比较如下图所示图2.2 基于叠加原理的打靶法算例图像图像中的直线表示边值问题的精确解，蓝色点表示利用算法1得到的数值解2.2 有限差分法有限差分方法是用于微分方程定解问题求解的最广泛的数值方法，其基本思想是用离散的、只含有有限个未知量的差分方程去近似代替连续变量的微分方程和定解条件，并把相应的差分方程的解作为微分方程定解问题的近似解。本节我们依然讨论边值问题
16、（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）可得一阶导
17、数的差分近似表达式为 （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分成等
18、份。设其节点分别为，。其中为步长，为内节点，称为边界节点。在每一个节点处，边值问题的方程（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 其他边值条件的有限差分方程在实际的应用中，我们可能会碰到
19、边值条件为（2.3）或（2.4）的边值问题，即第二或第三边值条件.对于第二或第三边值条件，在建立差分方程时，必须选用适当的差商公式代替边界上的和。例如，为了使近似解在边界和上具有二阶精度，可以选用本章第一节表2.2中列出的差商公式近似： （2.38）在（2.31）中，取近似那么我们可以得到第二边值条件的差分方程组: （2.39）第三边值条件的差分方程组： （2.40）其中包含个未知量。其系数矩阵都具有下列形态.其中，*代表不恒为零的元素。由线性代数的知识可以知道，矩阵可经过适当的初等行变换可化为三对角矩阵。2.2.4 有限差分方程的解法下面我们来研究下这些差分方程组的解法（以第一边值条件所得的
20、差分方程组（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,33.1 计算for i=2，3，N-2 3.2 解 （该存储的值）for i=2，3，N-13.3 解（为了节省空间用继续用存储（2.41）的解）for i=N-2，2，14 输出相应的值.例2.2
21、 用有限差分法解例题2.1中的边值问题：考虑有限差分法来求解这个问题。在算法2中，取，。计算结果见下表2.3及其图像见图2.3，值逼近，表示逐点绝对误差。表 2.3 有限差分法算例数据1.092600520-0.2877776121.187043128-0.4171168861.283336870-0.4549386471.381402046-0.4390546841.4811202621.481159417-0.3915488761.5823598951.582392461-0.3256506291.684989018-0.2494334951.7888817461.788898535-0.1678844821.893921099-0.084100540精确解与用有限差分法求得的数值解的比较如下图所示图 2.3 有限差分法算例图像图像中的直线表示边值问题的精确解，蓝色点表示由算法2得到的数值解.从该题的结果中我们可以看出，例2.2的数值解的精度与例2.1相比差很多。原因是在例2.1中所用的方法是建立在Runge-Kutta算法的基础上的，而Runge-Kutta算法具有的精度，而有限差分法的局部
