微分方程数值解及Matlab实现(1)

《微分方程数值解及Matlab实现(1)》由会员上传分享，免费在线阅读，更多相关内容在行业资料-天天文库。

1、铜陵职业技术学院学报2011年第3期微分方程数值解及Matlab实现王宝珍1,2（1.郑州大学，河南郑州450000；2.黄淮学院，河南驻马店463000）摘要：实际问题和科学研究中所遇到的微分方程往往很复杂，很多情况下不可能求出它的解析解，只需要数值解，分析了龙格-库塔法的基本原理，并给出了相应步骤、Matlab程序，针对目前比较前沿的脉冲微分方程与时滞微分方程的数值算法进行了设计，运用经典的四阶龙格库塔方法给出了脉冲微分方程和时滞微分方程以及脉冲时滞微分方程数值解法的计算步骤与相应的程序实现.

2、关键词：微分方程；数值解；Matlab实现；仿真实例中图分类号：O175文献标识码：A文章编号：1671-152X(2011)03-0095-03微分方程在经济、生物科学、化工等科技领域中得到了普进式”，即求解过程顺着节点排列的次序一步一步地向前推进.遍的应用，然而许多实际问题和科学研究中所遇到的微分方描述这类算法，只要给出从已知信息yn,yn-1,yn-2,…来计算yn+1的程往往很复杂，有时只需要数值解，因此对微分方程数值解求递推式，这类计算公式统称差分格式。解方法的理论整合及层次性的分析，并

3、结合数学软件Matlaby（x）-y（x）如果我们考察差商n+1n，根据微分中值定理，存把这些理论方法加以实现是有必要的。基于泰勒级数构造出h的龙格-库塔法，特别是四阶龙格—库塔法，其优点是精度高、在点ζ，x＜ζ＜x使得y（xn+1）-y（xn）=y（1ζ），从而利用所给方程nn+1程序简单、计算过程稳定，且易于调节步长。针对目前比较前h得y1=f，得：沿的脉冲微分方程与时滞微分方程的数值算法进行了设计，y（x）=y（x）+hf（ζ,y（ζ））运用经典的四阶龙哥库塔方法给出了脉冲微分方程和时滞微n

4、+1n其中的k*=f（ζ,y（ζ））称作区间[x,x]上的平均斜率。这样，只分方程以及多维微分方程数值解法的计算步骤与相应的程序nn+1实现，最后给出了仿真实例，证明了该算法是可行的与高效的。要对平均斜率k*提供一种算法，由上式便相应地导出一种计1.龙格—库塔法的基本原理算格式。科学技术中常常需要求解常微分方程的定解问题，这类按照这种观点考察欧拉格式，它简单地取点的斜率值k1=问题的最简单形式，是本章将要着重考察的一阶方程初值问（x,y）作为区间[x,x]上的平均斜率k*，精度自然很低。nnnn+

5、1题的数值解法。再考查改进的欧拉格式，它可改写成下列平均化形式：y'=f(x,y)"h$y=y+（k+k）!$n+1n12y(x)=y$200$$#）$k=f（x,y我们假定右函数适当光滑，譬如关于满足李普希兹$1nn$$（Lipshitz）条件，以保证上面常微分方程初值问题的解存在且$k=f（xy+hk）%2n+1,n1唯一。因此可以理解为，它用x与x两个点的斜率值k和k取nn+112虽然求解常微分方程有各种各样的解析方法，但解析方算术平均作为平均斜率k*，而x处的斜率值则利用已知信息n+1法只

6、能用来求解一些特殊类型的方程，求解从实际问题当中yn通过欧拉方法来预报。归结出来的微分方程主要靠数值解法。这个处理过程启示我们，如果设法在[x,x]内多预报几nn+1差分方法是一类重要的数值解法。这类方法是要寻找一个点的斜率值，然后将它们加权平均作为区间[x,x]上的平nn+1系列离散节点x1＜x2＜…＜…＜xn上的近似解y1,y2,…，yn…。相邻均斜率k*，则有可能构造出更高精度的计算格式，这就是龙两个节点的间距h=xn+1-xn称为步长。今后如不特别申明，总是格—库塔（Runge—Kutta

7、）方法的设计思想。假定为常数。为了使算法具有较高的精度，下面就作进一步讨论，随意初值问题的各种查分方法有个基本特点，它们都采取“步收稿日期：2011-05-24作者简介：王宝珍（1979-），女，河南南阳人，郑州大学数学系在职硕士研究生在读，黄淮学院数学系讲师，主要从事计算数学研究。·95·考察区间[x,x]内一点:2.龙格—库塔法模拟三维Lorenz方程nn+1x=x+ph,0＜p≤1，Lorenz方程是一个微型大气模型的简化形式，是麻省理n+pn我们希望用x，x两点的斜率值k1和k2加权平均得

8、到斜率工学院的气象学家E.Lorenz设计的这个模型来研究Rayleihnn+pBenard对流，这是热在流体（例如空气）中从较低的热介质（例k*，即令如地面）到较高的冷介质（例如上面的大气层）的运动。在这个yn+1=yn+h（[1-λ）k1+λk2]二维大气模型中形成了空气循环，这可用以下的三个方程的式中λ为待定系数，这里取k=f（x,y），问题在于该怎样预报x1nnn+p方程组来描述：处得斜率k2。≤x＇=-sx+sy≤≤我们先用欧拉方法提供y（x）的预报值y：≤＇n+pn+p