毕业论文---微分方程数值解算法分析与Matlab实现.doc

微分方程数值解算法分析与Matlab实现摘 要随着社会的发展，微分方程在经济、生物科学、化工等许多科技领域中得到了越来越广泛的应用，而所建立的微分方程模型中，只有极少数的微分方程可以得到解析解，为了了解微分方程的解的性质，不得不进行微分方程数值解的探讨与研究，所以对微分方程数值解的研究有重要的应用价值.本文首先讨论了常微分方程数值解的几种常见的算法(如欧拉法、改进欧拉法、龙格库塔法等)，分析了算法的具体设计思想，并给出了相应步骤、Matlab程序，及其对算法的误差做了相应的分析，通过实例验证了算法的可行性与有效性.其次对于部分算法的收敛性与稳定性从理论上给出了分析，证明上述算法是收敛的与稳定的.最后针对目前比较前沿的脉冲微分方程与时滞微分方程的数值算法进行了设计，运用经典的四阶龙格库塔方法给出了脉冲微分方程和时滞微分方程以及脉冲时滞微分方程数值解法的计算步骤与相应的程序实现，最后给出了测试实例，证明了该算法是可行的与高效的.关键词：微分方程；数值解；Matlab实现；仿真实例 AbstractWith the development of society, differential equations have been used widely in the economic, biological sciences, chemical and many other fields of technology .But in the differential equation model which was established, only very few analytical solutions of differential equations can be obtained, In order to understand the nature of the solution of differential equations, we have to explore and research differential equations ,so the research on the numerical solution of differential equations have important applications.Firstly，the numerical solution of ODEs of several common algorithms (such as Euler, improved Euler, Runge-Kutta, etc) is discussed, the theoretical algorithms are analyzed in detail, the relevant steps, Matlab programs, and the corresponding analysis of the errors of the theoretical algorithms are given. The feasibility and effectiveness of the algorithms are proved by solving examples. Secondly, the convergence and stability analysis of parts of algorithms are given in theory, the above algorithms are convergence and stable. Finally, the numerical algorithms of the relatively cutting-edge Impulsive differential equations and delay differential equations are designed, the calculation steps, the corresponding program implementation of the delay differential equations ,the impulsive differential equations and the numerical solution of impulsive and delay differential equations are given by using classical Runge-Kutta methods, it shows that the algorithms are feasible and efficient by giving some examples of numerical simulation in the last.Key words：Differential equation; Numerical solution; Matlab simulation; ExamplesIIII目 录摘 要IAbstractII1 引 言12 微分方程数值解算法22.1 欧拉方法22.2 向后欧拉格式32.3 梯形格式52.4 改进欧拉格式62.5 截断误差分析72.6 龙格-库塔法92.7 亚当姆斯方法133 差分方法的收敛性和稳定性163.1 关于差分方程163.2 差分方法的相容性、收敛性163.2.1 单步法的相容性163.2.2 单步法的收敛性173.3 差分方程的稳定性183.3.1 差分方法的绝对稳定性183.3.2 数值算法的稳定性分析204 脉冲时滞微分方程的数值解法234.1 Runge-Kutta法求脉冲时滞微分方程的算法分析234.2 求解脉冲时滞微分方程的算法步骤244.3 龙格库塔法求解脉冲和时滞微分方程的仿真实例254.4 龙格库塔法求解脉冲时滞微分方程的仿真实例265 结束语28参考文献29致 谢301 引 言微分方程在解决经济、生物科学、化工等科技领域中得到了普遍的应用，许多实际问题的数学模型是微分方程或者微分方程组的定解问题,如物体运动、电路振荡、化学反应及生物群体的变化等.实际问题和科学研究中所遇到的微分方程往往很复杂，很多情况下不可能求出它的解析解，有时候即使求出它的解析解，也会由于解得表达式过于繁复而不实用，并且许多实际问题不需要求解微分方程的解析解，只需要数值解，因此利用数值解法求解实际问题具有十分重要的理论意义及使用价值. 现如今许多文献对某些特殊微分方程的解析算法进行分析，但缺乏对微分方程数值解的常见算法的系统性的理论分析与整合，使得此方面的知识得不到系统的认识.然而对微分方程数值解求解方法的理论整合及层次性的分析，并结合数学软件Matlab把这些理论方法加以实现是有必要的，既能在理论上将微分方程数值解的求法得到了系统性的研究，又能将这些理论知识更方便地运用到实际问题中.本文主要探讨的是微分方程数值解的算法和Matlab实现的具体问题. 内容主要是介绍求解微分方程的初值问题的差分方法，并通过对数值解算法的严密分析以及对实际应用例题的解法综述，来表明算法的精确性和算法程序化的可行性由于微分方程有多种类型，其特征也不能逐一而论，故其算法也具体多样，本文对微分方程常见的算法如欧拉方法、改进的欧拉方法、龙格-库塔方法及亚当姆斯方法等作出系统的分析和总结，本文首先介绍的是欧拉方法，其中包括了欧拉格式、向后欧拉格式、梯形格式、改进的欧拉格式等.欧拉格式形式简单，便于理解，但是精度不高，所以提出了向后欧拉格式和梯形格式，以提高欧拉公式的精度.其次，本文介绍了基于泰勒级数构造出的龙格-库塔法，特别是四阶龙格-库塔法，其优点是精度高、程序简单、计算过程稳定，且易于调节步长.龙格库塔方法是一类重要算法，但这类算法在每一步都需要先预报几个点上的斜率值，计算量比较大.基于这一缺陷，本文提出了亚当姆斯法，并对亚当姆斯法的预报-校正系统作出算法分析.最后，对于部分算法的收敛性与稳定性从理论上给出了分析，证明上述算法是收敛的与稳定的.最后针对目前比较前沿的脉冲微分方程与时滞微分方程的数值算法进行了设计，运用经典的四阶龙哥库塔方法给出了脉冲微分方程和时滞微分方程以及多维微分方程数值解法的计算步骤与相应的程序实现，最后给出了测试实例，证明了该算法是可行的与高效的.2 微分方程数值解算法科学技术中常常需要求解常微分方程的定解问题.这类问题的最简单形式，是本章将要着重考察的一阶方程初值问题的数值解法. (2-1)我们假定右函数适当光滑，譬如关于满足李普希兹(Lipshitz)条件12，以保证上面常微分方程初值问题的解存在且唯一.虽然求解常微分方程有各种各样的解析方法，但解析方法只能用来求解一些特殊类型的方程，求解从实际问题当中归结出来的微分方程主要靠数值解法.差分方法是一类重要的数值解法.这类方法是要寻找一系列离散节点上的近似解.相邻两个节点的间距称为步长.今后如不特别申明，总是假定为常数.初值问题的各种查分方法有个基本特点，它们都采取“步进式”，即求解过程顺着节点排列的次序一步一步地向前推进.描述这类算法，只要给出从已知信息来计算的递推式，这类计算公式统称差分格式12.2.1 欧拉方法2.1.1 欧拉格式的基本原理方程中含有导数项，这是微分方程的本质特征，也正是它难以求解的症结所在数值解法的第一步就是设法消除其导数项，这项手续称离散化由于差分是微分的近似计算，实现离散化的基本途径是用差商代替导数，譬如，若在点处列出方程并用差商替代其中的导数项，结果有设用的近似值代入上式的右端，记所得结果为，这样导出的计算公式 (2-2)称为欧拉公式2.1.2 欧拉格式的Matlab实现function E=Euler_1(fun,x0,y0,xN,N)%欧拉公式，其中：fun为一阶微分方程的函数；x0,y0为初始条件； xN为取值范围的一个端点；h为区间步长；N为区间的个数；x为Xn构成的向量； y为Yn构成的向量x=zeros(1,N+1);y=zeros(1,N+1);x(1)=x0;y(1)=y0;h=(xN-x0)/N;for n=1:Nx(n+1)=x(n)+h;y(n+1)=y(n)+h*feval(fun,x(n),y(n);endT=x,y2.1.3 欧拉格式实例分析例2.1 求解初值问题解 是用欧拉方法，建立一个M-文件，其内容为图 2.1 欧拉法运行结果示意图function z=f(x,y)z=y-2*x/y;取迭代步骤为10.在Matlab运行窗口中输入： E=Euler_1(f,0,1,1,10)输出结果如图(2.1)所示：2.2向后欧拉格式2.2.1 向后欧拉格式的基本原理如果在微分方程离散化时，用向后差商代替导数，即则得到计算公式 (2-3)用这组公式求问题式(2-1)的数值解称为向后欧拉公式.向后欧拉公式与欧拉公式形式上相似，但实际计算时却复杂得多.欧拉公式是显式的，可直接求解.向后欧拉公式(2-3)第一式的右端含有，因此是隐式公式，则称向后欧拉公式为单步隐式公式，一般要用迭代法求解，迭代公式通式为(2-4)2.2.2 向后欧拉格式的Matlab实现function E2=Euler_2(fun,x0,y0,xN,N)%欧拉公式，其中：fun为一阶微分方程的函数；x0,y0为初始条件； xN为取值范围的一个端点；h为区间步长；N为区间的个数；x为Xn构成的向量； y为Yn构成的向量x=zeros(1,N+1);y=zeros(1,N+1);x(1)=x0;y(1)=y0;h=(xN-x0)/N;for n=1:N%用迭代法求 x(n+1)=x(n)+h; z0=y(n)+h*feval(fun,x(n),y(n);for k=1:3 z1=y(n)+h*feval(fun,x(n+1),z0); if abs(z1-z0) Euler_2(f,0,1,1,10)输出的结果如表(2.1)其结果如图(2.2).x00.10000.20000.30000.40000.50000.60000.70000.80000.90001.0000y1.00000.90910.82640.75120.68280.62070.56420.51290.46620.42380.3852程序计算得近似值为,而精确值为3误差较大.梯形公式(下一算法)与向后欧拉公式相比较，梯形公式的算法精确度较高.2.3 梯形格式为了提高精度，我们改用梯形方法计算积分项代入式(2-3)，有 设将式中，的分别用，替代，作为离散化的结果导出下列计算格式 (2-5)上述差分格式称为梯形公式.梯形公式是隐式格式，一般需用迭代法求解，迭代公式2为：(2-6)2.3.1 梯形格式的Matlab实现function E3=Euler_3(fun,x0,y0,xN,N)%梯形公式，其中：fun为一阶微分方程的函数；x0,y0为初始条件；xN为取值范围的一个端点；h为区间步长；N为区间的个数；x为Xn构成的向量；y为Yn构成的向量x=zeros(1,N+1);y=zeros(1,N+1);x(1)=x0;y(1)=y0;h=(xN-x0)/N;for n=1:Nx(n+1)=x(n)+h;z0=y(n)+h*feval(fun,x(n),y(n);for k=1:3 z1=y(n)+h/2*(feval(fun,x(n),y(n)+feval(fun,x(n+1),z0); if abs(z1-z0) Euler_3(f,0,1,1,10)输出结果如图(2.3)所示：程序计算得近似值为，而精确值为3误差较小.梯形公式与欧拉公式相比较，梯形公式的算法精度有所提高. 2.4 改进欧拉格式 综合使用欧拉方法和梯形方法，即先用欧拉方法求得一个初步的近似值，记为，称之为预报值；预报值的精度不高，我们用它替代式子右端再直接计算，得到校正值，这样建立的预报-校正系统 (2-7)称作改进的欧拉公式.这是一种一步显式格式，它可表为如下嵌套形式或表为下列平均化形式1 2.4.1 改进的欧拉格式的Matlab实现function E4=Euler_4(fun,x0,y0,xN,N)%改进的欧拉公式，其中：fun为一阶微分方程的函数；x0,y0为初始条件；xN为取值范围的一个端点；h为区间步长；N为区间的个数；x为Xn构成的向量；y为Yn构成的向量x=zeros(1,N+1);y=zeros(1,N+1);x(1)=x0;y(1)=y0;h=(xN-x0)/N;for n=1:Nx(n+1)=x(n)+h;z0=y(n)+h*feval(fun,x(n),y(n);y(n+1)=y(n)+h/2*(feval(fun,x(n),y(n)+feval(fun,x(n+1),z0);endT=x,y2.4.2 改进的欧拉格式实例分析例2.4 用改进的欧拉方法求解初边值问题图 2.4 改进欧拉法运行结果示图的解在x=0.1处得近似值，迭代次数N=5.解 建立一个M-文件，其内容为function z=f(x,y)z=x+y;取迭代次数为5.在Matlab运行窗口中输入： Euler_4(f,0,1,0.1,5)输出结果如图(2.4)所示：2.5 截断误差分析欧拉公式在计算时只用到了前一步的值。另外，若已知，欧拉公式给出了，的显式依赖关系，将代入(2-2)右端可直接得到。我们称欧拉公式的单步显示公式。单步显式公式的一般形式为(2-8)称为增量函数。欧拉公式(2-2)的增量函数为。一般来说微分方程问题的精确解不满足式(2-8)，即一般地有2.5.1 局部截断误差定义2.1 称为单步显式公式(2-8)在处得局部截断误差。一个求解公式的局部截断误差刻画了其逼近微分方程的准确程度。根据上述定义可直接求得欧拉公式(2-2)的局部截断误差还可以求得梯形公式的隐形格式(2-5)的截断误差为2.5.2 整体截断误差用某种数值方法(例如Euler方法，改进Euler方法)求得的数值解，一般来说与步长有关，为了反映出这种关系，我们将其记为，.求数值解的目的就是用作为的近似值.人们自然要问在每一节点处近似值与精确解的差是否很小？定义2.2 设为微分方程初值问题的解在节点处的值，为用某种数值方法求得的近似解.称 (2-6)为该方法的整体截断误差.如果则称该方法是收敛的.整体截断误差为所有节点上误差的最大值.他和局部截断误差是有紧密联系的.在一定的条件下，如果局部截断误差是，则整体截断误差是.而分析局部截断误差是比较容易的，所以我们直接根据局部截断误差来刻画求解公式的精度，为此给出下面的定义.定义2.3 如果一个求解公式的局部截断误差为，则称该求解公式是阶的，或具有阶精度.根据此定义，Euler公式和向后Euler公式是一阶的，梯形公式和改进Euler公式是二阶的.本章所介绍的Euler公式、向后Euler公式和梯形公式是由数值积分得到的，改进Euler公式是由预测校正方法得到的.事实上，这些公式以及具有更高精度的求解公式还可用其它方法推得.2.6 龙格-库塔法考察差商，根据微分中值定理，存在点，使得从而利用所给方程得， (2-7)其中的称作区间上的平均斜率.这样，只要对平均斜率提供一种算法，由上式便相应地导出一种计算格式.按照这种观点考察欧拉格式，它简单地取点的斜率值作为区间上的平均斜率，精度自然很低.再考查改进的欧拉格式，它可改写成下列平均化形式12：因此可以理解为，它用与两个点的斜率值和取算术平均作为平均斜率，而处的斜率值则利用已知信息通过欧拉方法来预报.这个处理过程启示我们，如果设法在内多预报几个点的斜率值，然后将它们加权平均作为区间上的平均斜率，则有可能构造出更高精度的计算格式，这就是龙格-库塔(Runge-Kutta)方法的设计思想.2.6.1 龙格-库塔法的基本原理为了使算法具有较高的精度，下面就作进一步讨论，随意考察区间内一点， 我们希望用，两点的斜率值和加权平均得到斜率，即令式中为待定系数，这里取，问题在于该怎样预报处得斜率？我们先用欧拉方法提供的预报值：然后用通过计算产生斜率值这样设计出的计算格式具体形式 (2-8) 其中含有两个待定参数，我们希望适当选取这些参数的值，使得格式(2-8)具有较高的精度.仍然假定，分别将和进行泰勒展开，有代入式(4-1)得和二阶泰勒展开式比较系数即可发现，欲使格式(2-7)的截断误差为，只要成立.满足这一条件的一族格式(2-7)统称二阶龙格-库塔格式.为了进一步提高精度，用类似二阶龙格-库塔的构造方法可以确定三阶和四阶龙格-库塔法的参数，来构造出三阶和四阶的龙格-库塔法.四阶龙格-库塔法也不只一个，下面给出常用的四阶经典的龙格-库塔公式: (2-9)2.6.2 龙格-库塔法的Matlab实现function R=rk4(f,a,b,ya,N)%y=f(x,y);a,b左右端点；N为迭代步长；h为步长；ya为初值；h=(b-a)/N;T=zeros(1,N+1);Y=zeros(1,N+1);t=a:h:b;Y(1)=ya;for j=1:Nk1=h*feval(f,T(j),Y(j);k2=h*feval(f,T(j)+h/2,Y(j)+k1/2);k3= h*feval(f,T(j)+h/2,Y(j)+k2/2);k4= h*feval(f,T(j)+h,Y(j)+k3);Y(j+1)=Y(j)+(k1+2*k2+2*k3+k4)/6;endR=T Y2.6.3 龙格-库塔法实例分析例2.5 用程序求解初值问题：图 2.5 龙格-库塔法运行结果示图解 建立一个M-文件，其内容为function z=f(x,y)z=y-2*x/y;在Matlab运行窗口中输入：rk4(f,0,1,1,10)输出结果如表(2.5)所示：并输出结果如图（2.5）T00.10000.20000.30000.40000.50000.60000.70000.80000.90001.0000Y1.00001.09541.18321.26491.34161.41421.48321.54921.61251.67331.7321表 2.5 龙格-库塔法运行结果2.6.4 龙格-库塔法模拟三维Lorenz方程Lorenz方程是一个微型大气模型的简化形式，是麻省理工学院的气象学家E.Lorenz设计的这个模型来研究Rayleih Benard对流，这是热在流体(例如空气)中从较低的热介质(例如地面)到较高的冷介质(例如上面的大气层)的运动.在这个二维大气模型中形成了空气循环，这可用以下的三个方程的方程组来描述11： (2-10)图 2.6 洛伦茨方程模拟仿真结果示意图变量表示顺时针循环的速度，表示上升和下降空气柱之间的温度差，而度量了来自垂直方向的严格线性温度剖面的偏差.Prandtl数s、Reynolds数及时方程组的参数.在大多数情况下，取参数为.应用微分方程数值解法的龙格库塔算法程序代码如下：function z=ydot(t,y)%Lorenz equationsS=10;r=28;b=8/3;z(1)=-s*y(1)+s*y(2);z(2)=-y(1)*y(3)+r*y(1)-y(2)z(3)=y(1)*y(2)-b*y(3)计算机模拟仿真结果如图(2.6)：2.7 亚当姆斯方法2.7.1 线性多步法 在逐步推进的求解过程中，计算之前事实上已经求出了一系列的近似值.如果充分利用前面多步的信息来预测，则可期望获得较高的精度，这就是构造多步法的基本思想.线性步法的一般公式为 (2-11)其中均为与无关的常数，.当时为显格式；当时为隐格式.特别当时为Euler公式；当时为梯形公式.下面我们将介绍基于Taylor展开的方法构造Adams法.2.7.2 亚当姆斯公式的基本原理 考察式子(2-7).设用两点的斜率值加权平均作为区间上的平均斜率，这样设计出的计算格式具有形式 (2-12)我们希望适当选取参数，使以上格式具有二阶龙格库塔方法的精度相当. 假设，将式子泰勒展开得由此可见，为使计算格式具有二阶精度，须取，这样导出的计算格式称作二阶亚当姆斯格式.类似地可以导出三阶亚当姆斯格式和四阶亚当姆斯格式这里和下面均记为上述亚当姆斯格式都是显式的，但由于用节点的斜率值来预报区间上的平均斜率是个外推过程，效果不够理想为了进一步改善精度，我们增加节点的斜率值来得出上的平均斜率譬如，考察形如 (2-13)的隐式格式，设，将上式右端泰勒展开有可见欲使格式具有二阶精度，需令，这样构造出的二阶隐式亚当姆斯格式其实就是梯形格式，类似地可以导出三阶隐式亚当姆斯格式和四阶亚当姆斯格式.Error! Reference source not found.2.7.3 亚当姆斯预报-校正系统 仿照改进欧拉公式的构造方法，将显式和隐式两种亚当姆斯格式相匹配，构成下列亚当姆斯预报-校正系统1 (2-14)亚当姆斯预报-校正系统是四步法，它在计算时不但要用到前一步的信息，而且要用到更前面三步的信息，因此它不能自行启动，在实际计算时，可借助于某种单步法，譬如用四阶龙格-库塔格式(2-9)为预报-校正系统(2-14)提供开始值，.2.7.4 亚当姆斯预报-校正系统的Matlab实现求预报-校正系统的微分方程问题的算法具体描述如下：1) 设定的初始值,用经典格式(2-9)提供，的开始值；2) 设，将n的值设定为4.3) 然后作如下赋值：；，并输出，；4) 判断是否等于,如果相等，则结束；如果不相等，则令；，再重复运行第三步，直到满足；5) 输出结果并结束.2.7.5 亚当姆斯法的实例分析例2.6 用亚当姆斯预报-校正系统(2-13)求解初值问题解 取步长，用四阶龙格-库塔格式(2-9)计算开始值，然后套用预报-校正系统(2-13)，计算结果见下表(2.6)：0.01.00001.00000.11.09541.09540.21.18321.18320.31.26491.26490.41.34151.34161.34160.51.41411.41421.41420.61.48321.48321.48320.71.54911.54921.54920.81.61241.61241.61250.91.67331.67331.67331.01.73201.73201.7321表 2.6 亚当姆斯法运行结果3 差分方法的收敛性和稳定性3.1 关于差分方程为了便于讨论差分方法的收敛性以及稳定性，本文给出了差分方程的一般理论.具有常系数的线性差分方程一般形式为 (3-1)为差分方程的起始下标.不失一般性，可设，当时，称(3-1)为阶差分方程，而 (3-2)称为(5-1)的齐次差分方程.若给定个初始值，(3-1)唯一地决定了差分方程解.对齐次方程，其解也可由个初始值唯一决定.3.2 差分方法的相容性、收敛性3.2.1 单步法的相容性常微分方程初始问题 (3-3)的单步法 (3-4)称(3-4)与常微分方程(3-3)是相容的，如果微分方程的理论解在极限情况下满足差分方程，即，由于因此，相容性条件为 (3-5)相容性可理解为，差分方程(3-4)确是微分方程(3-3)的近似，而不是其它方程的近似.容易验证，我们导出的单步方法，诸如Euler法，Euler向后公式，Runge-Kutta法均为相容方法.事实上，对Euler法而向后Euler法故也有.对Runge-Kutta法为了避免符号理解上的困难，不妨记其中.显然有故有，但在Runge-Kutta方法推导过程中，系数满足的第一个方程就是因此，即Runge-Kutta法是相容格式.3.2.2 单步法的收敛性关于单步法的收敛性，本论文前面已有所讨论.由导出的单步法当在上满足Lipschitz条件是，差分方程数值解误差满足 (3-6)因此，若，则单步方法是收敛的，以此得到如下定理6定理 3.1 若微分方程(3-3)由定义的单步方法对，若存在，当时，在上连续并满足Lipschitz条件则单步法收敛的充分必要条件为格式是相容的.证明 充分性：若满足定理条件时，单步法误差满足当单步法相容，则单步法是微分方程近似，故有，其中，因此，单步法收敛.必要性：若，则单步法收敛时，只能收敛于微分方程 (3-7)的解.但，(3-3)与(3-5)决定的是两个不同的曲线，故单步法不可能收敛于(3-3)的解.3.3 差分方程的稳定性对差分方程进行数值求解时，由于初始值可能带有误差，以及每一步计算时的舍入误差，因此，数值结果与差分方程理论解之间存在误差.而且这一步产生的计算误差将带到以后的计算中，这类误差是否能得到控制，它与差分方程本身有关，能控制误差的差分方程为稳定格式，否则为不稳定格式.这就需要对差分方程作稳定性分析，差分方程的渐进稳定性与差分方程的收敛性等价，本文重点分析了差分方程的绝对稳定性.3.3.1 差分方法的绝对稳定性定义3.2 差分方程 (3-7)称为绝对稳定的，如果差分方程作用到微分方程时，对任意初始值，总存在左复平面的一个区域，当取自这个区域时，差分方程的解趋于0.而这个区域称为稳定区域.当稳定区域为整个左半复平面时，称差分方程为无条件稳定的.定理3.2 差分方程(3-7)绝对稳定的充分必要条件为，存在左半复平面上的一个区域时，齐次差分方程6 (3-8)的特征多项式的根按模小于1.证明 差分方程(3-7)作用到时，成为.上式为齐次差分方程，其特征多项式为 (3-9)对任给的一组初始值，差分解趋于零的充分必要条件为(3-8)的特征方程(3-9)的特征根按模小于1，而(3-9)的根为的函数，因此，若存在左半复平面上区域，取自这个区域时，(3-9)的特征根按模小于1，则按绝对稳定定义，差分方程(3-8)时绝对稳定的.相比于渐进稳定性，绝对问法定性更有现实意义，不满足绝对稳定的差分方程，由于不存在稳定区域，因此找不到适当的，使差分方程的解对初始误差有所抑制，这样的格式就不能使用，而绝对稳定的格式中存在，当时，差分解受初始误差的影响随着计算一步步推进而逐渐消失.当然，的大小也决定于的界()，的界相当于.下面将对各类差分方法详细讨论其绝对稳定性.3.3.2 数值算法的稳定性分析讨论Euler公式yx-10-2图 3.1 Euler法的稳定区域的绝对稳定性.Euler公式作用于时，成为其特征方程的特征根为当落在以(-1,0)为圆心的单位圆时，因此Euler法为绝对稳定的，稳定区域如图(3.1)所示.讨论向后Euler公式的绝对稳定性.向后Euler公式作用于后有其特征方程为的特征根为当时，恒有，故向后Euler公式为无条件绝对稳定格式.以三阶龙格库塔法为例的龙格-库塔法的绝对稳定性.格式作用到，有故也有其特征方程为而特征根为对，存在，当时，特征根满足因此三阶龙格库塔法为绝对稳定格式.事实上，我们可以证明，阶的龙格库塔法作用于时，有故有特征多项式于是，对，存在，当时，有图 3.2 R-K法稳定区域亦即，龙格库塔法均为绝对稳定格式.而且同阶的不同格式有相同的稳定区域，但不同阶的格式，它们的稳定区域也不同.图(3.2)给出了各阶龙格库塔法的稳定区域.再来讨论Adams公式的绝对稳定性步Adams显格式作用于，有它的特征方程为当时，有特征根由于特征方程根连续依赖于，因此当充分小时，仍有.记并代入特征方程，则有两端展开后比较的同阶系数，特别由的一次系数，可得到故有 由于为上数值积分的积分系数和，因而，于是当充分小是，Adams显格式绝对稳定.类似地，也可证明Adams隐格式的绝对稳定性.特别，时的向后欧拉公式及时的梯形公式是无条件绝对稳定的.4 脉冲时滞微分方程的数值解法4.1 Runge-Kutta法求脉冲时滞微分方程的算法分析考虑如下脉冲时滞微分方程初值问题10： (4-1)求初值问题(4-1)的数值解，就是寻求准确解在一系列离散节点上的近似值，其中：脉冲时刻满足，当时，；是脉冲算子；是正整数的集合；是时滞；是初始值函数.显然，只有当初值问题(4-1)的解存在且唯一时，使用数值解法才有意义.当是连续且满足Lipchitz条件时，Lipchitz常数为，即对任意恒有时初值问题的解存在并且唯一.下面将脉冲时滞微分方程初值问题用四级四阶龙格库塔法表示出来，其截断误差为，每步都要计算4次的值，其标准公式即古典Runge-Kutta法为： (4-2)其中表示计算过程中选取的步长，通常采用等距节点.为了得到系统的精确解，一般步长选的尽可能小.若，则将区间由右向左以步长h进行分割.设则代表网域，其中元素称为节点.(4-2)式中的表示点处得斜率；表示利用求得的点处得斜率；表示利用求得的点处的斜率；表示利用求得的点处的斜率9.(4-2)式是常微分方程的数值解的龙格库塔法，而所谓时滞，就是一种时间差，是时间延迟，是控制系统普遍存在的一种现象，与始点和终点都有关，即某一时刻系统(4-1)的数值解与系统(4-1)在时刻及时刻的数值解都有关.将(4-2)式应用于方程(4-1)，得 (4-3)即当计算的值时，需要用到的值()，由于的值也在网格上，因此将其代替时滞项的值进行迭代运算.所谓脉冲，就是在短时间内突变，随后又迅速返回其初始值的物理量，他可以是周期性重复的，也可以是非周期性的或单次的.即在某一个脉冲时刻，将系统的解修改为所给定的脉冲值，而在该脉冲时刻和下一个脉冲时刻之间系统被重新积分.4.2 求解脉冲时滞微分方程的算法步骤求脉冲时滞微分方程初值问题的算法具体描述如下：1) 设定的初始值，并设定存储空间为；2) 设，应用(4-2)式分步计算区间的历史函数值；3) 设脉冲时刻为，则在每一个脉冲时刻，将该时刻的函数值修改为脉冲值，即将该时刻对应的存储单元的值修改为，再应用(4-3)式分步计算区间的系统(4-1)的数值解.求脉冲微分方程初值问题的算法具体描述如下：1) 设定的初始值，并设定存储空间为；2) 设，设脉冲时刻为，则在每一个脉冲时刻，将该时刻的函数值修改为脉冲值，即将该时刻对应的存储单元的值修改为，再应用(4-2)式分步计算区间的脉冲微分方程数值解.求时滞微分方程数值解的相应算法描述如下：1) 设定的初始值，并设定存储空间为；2) 设，应用(4-2)式分步计算区间的历史函数值；3) 在区域应用(4-3)式分步计算时滞微分方程的数值解.以上算法均采用Matlab编程实现.4.3 龙格库塔法求解脉冲和时滞微分方程的仿真实例例4.1 求如下脉冲微分方程的初值问题的数值解: (4-4)应用上述算法计算所得的数值解仿真结果如图(4.1)所示.图4.3 脉冲时滞微分方程仿真结果示意图图4.1 脉冲微分方程仿真结果示意图例4.2 对如下时滞Hopfield神经网络：(4-5)其中： ；，.应用上述算法得到其仿真结果如图(4.2)所示：图4.2 时滞微分方程仿真结果示意图4.4 龙格库塔法求解脉冲时滞微分方程的仿真实例例4.3 求如下脉冲时滞方程初值问题的数值解： (4-6)应用上述算法计所得的数值解仿真结果如图(4.3)所示例4.4 对如下脉冲时滞Hopfield神经网络：(4-7)其中： .应用上述算法得到其仿真结果如图(4.4)所示.图4.4 脉冲时滞Hopfield神经网络真结果示意图. 将龙格库塔法应用于求解脉冲微分方程和时滞微分方程，脉冲时滞微分方程以及多维微分方程的数值解，并运用Matlab编程实现，仿真结果表明该方法是高效的.由于脉冲时滞微分方程在计算机的辅助设计、电路分析、气象分析、力学系统、化学反应模拟及自动控制系统的实时仿真等科学与工程应用领域中有广泛的应用，因此其数值仿真方法对这些领域的实际科研工作有一定的参考价值. 5 结束语本文介绍求解微分方程的初值问题的方法.首先介绍的是欧拉方法，其中包括了欧拉格式、向后欧拉格式、梯形格式、改进的欧拉格式等.欧拉格式形式简单，便于理解，但是精度不高，所以有向后欧拉格式和梯形格式，以提高欧拉公式的精度.在实际应用中，常用方法是基于泰勒级数展开的构造方法和基于数值积分的构造方法.基于泰勒展开的方法更灵活，也具有一般性.另外，泰勒展开还有一个优点就是在构造公式的同时，可以得到关于截断误差的估计.而直接用泰勒展开导出的泰勒级数发不便于实际应用，所以本文重点介绍了基于泰勒级数构造出的龙格-库塔法，特别是四阶龙格-库塔法，是计算机上常用算法.四阶龙格-库塔法优点是精度高、程序简单、计算过程稳定，且易于调节步长.本文中给出了相应的Matlab程序，并采用四阶龙格-库塔法采用