常微分方程组初值问题数值解的实现和算法分析.doc

课程设计说明书（论文） 第 17 页 常微分方程组初值问题数值解的实现和算法分析摘 要本次课程设计主要内容是用改进Euler方法和四阶Runge-Kutta方法解决常微分方程组初值问题的数值解法，通过分析给定题目使用Matlab编写程序计算结果并绘图然后区别两种方法的使用范围。最后对计算结果进行分析，得到结论。关键词：改进Euler，Runge-Kutta，初值问题目 录1前言12题目叙述13解题思路13.1一阶常微分方程的初值问题13.2一阶常微分方程组的初值问题23.2.1用Runge-Kutta方法计算解决一阶微分方程组初值问题的基本思路23.2.2用改进Euler方法计算解决一阶微分方程组初值问题的基本思路44用matlab语言编程解决相关问题44.1四阶Runge-kutta方法的Matlab编程实现44.2 Euler改进方法Matlab编程实现55 编程解决65.1 输入计算题目65.2用Runge-Kutta方法的Matlab编程解法65.3用改进Euler方法的Matlab编程解法76计算结果86.1用四阶Runge-Kutta方法的Matlab编程解法的结果以及与精确解的比较96.2用改进Euler方法的Matlab编程解法的结果以及与精确解的比较97.结果分析10致谢11参考文献12附录13翻译171前言常微分方程是解决工程实例的常用的工具，建立微分方程只是解决问题的第一步，通常需要求出方程的解来说明实际现象，并加以检验。如果能得到解析形式的解固然是便于分析和应用的，但是我们知道，虽然求解常微分方程有各种各样的解析方法，但解析方法只能用来求解一些特殊类型的方程，实际问题终归结出来的微分方程主要靠数值解法1。数值解法就是一个十分重要的手段，而Euler方法以及Runge-Kutta方法又是数值解法最基础最常用的方法。通过对两种方法计算结果的对比，从而找到一个最为合适的计算方法。2题目叙述对以下的常微分方程组初值问题(1) 用四阶RungeKutta方法求解，根据计算结果画出解的图形(2) 用改进Euler方法求解，根据计算结果画出解的图形。3解题思路3.1一阶常微分方程的初值问题一阶常微分方程的初值问题的一般形式为： 我们知道，只要函数f (x, y)适当光滑譬如关于y 满足利普希茨(Lipschitz)条件，即如果存在实数，使得理论上就可以保证初值问题的解存在并且唯一。所谓数值解法， 就是求上述问题在一系列离散点的近似值。两个相邻节点的间距称为步长13.2一阶常微分方程组的初值问题解一阶微分方程组时类似于单个方程的数值解法，只要把和理解为向量，那么就有：一阶微分方程组初值问题的形式为： （1）式中引入向量符号： 则（1）可写为：(2)3.2.1用Runge-Kutta方法计算解决一阶微分方程组初值问题的基本思路(2)式形式上与常数微分方程初值问题是一样的，只要注意向量函数运算及其表示，就可以用初值问题的求解格式得到常微分方程组初值问题（2）的求解格式，由初值问题的经典Runge-kutta公式可得一阶常微分方程组初值问题（2）的Runge-kutta公式：注意上式是向量形式，其对应的分量形式为：微分方程理论告诉我们，高阶微分方程可转化为一阶微分方程组来研究，因此可以用一阶微分方程组初值问题揭发来解高阶微分方程初值问题。高阶微分方程初值问题的形式为： (3)令则（2）化为了一阶微分方程组初值问题：Runge-kutta方法巧妙利用函数在一些点上的函数值的线性组合，获得了高阶的数值解法，它避开了要获得高阶方法须对求高阶导数的不便，是离散化方法中Tayl情况，其中在准确性的工作量的综合效果看，经典的Runge-kutta方法是首选or展开法的一个应用。Runge-kutta方法主要用于定步长的。Runge-kutta方法也常用于对多步法提供初值。3.2.2用改进Euler方法计算解决一阶微分方程组初值问题的基本思路改进Euler方法需要用Euler方法求出一个预测值然后再用梯形公式校正一次得到，即所求结果的迭代格式。 （4）为了方便编程可将（4）式改变为如下格式4用matlab语言编程解决相关问题4.1四阶Runge-kutta方法的Matlab编程实现function T=Runge_Kutta(f,x0,y0,h,n)% T=Runge_Kutta(f,x0,y0,h,n)% f 待解方程(组)% x0 初试自变量值% y0 初试函数值% h 步长% n 步数if nargin5 n=100;endr=size(y0); r=r(1);s=size(x0); s=s(1);r=r+s;T=zeros(r,n+1);T(:,1)=y0;x0;for t=2:n+1k1=feval(f,T(1:r-s,t-1);x0);k2=feval(f,k1*(h/2)+T(1:r-s,t-1);x0+h/2);k3=feval(f,k2*(h/2)+T(1:r-s,t-1);x0+h/2);k4=feval(f,k3*h+T(1:r-s,t-1);x0+h);x0=x0+h;T(:,t)=T(1:r-s,t-1)+(k1+k2*2+k3*2+k4)*(h/6);x0;End 24.2 Euler改进方法Matlab编程实现用Euler改进方法编写Matlab程序为 ：function x,y=eulerpro(fun,x0,xfinal,y0,n);h=(xfinal-x0)/n;x(1)=x0;y(1)=y0;for i=1:nx(i+1)=x(i)+h;y1=y(i)+h*feval(fun,x(i),y(i);y2=y(i)+h*feval(fun,x(i+1),y1);y(i+1)=(y1+y2)/2;end5 编程解决5.1 输入计算题目function dY=dydx(X,Y)dY(1)=-0.013*Y(1)-1000*Y(1)*Y(2);dY(2)=-2500*Y(2)*Y(3);dY(3)=0.013*Y(1)-1000*Y(1)*Y(2)-2500*Y(2)*Y(3) dY=dY(1);dY(2);dY(3);5.2用Runge-Kutta方法的Matlab编程解法function k,X,Y,wucha,P=RK4z(dydx,a,b,CT,h)n=fix(b-a)/h);X=zeros(n+1,1); Y=zeros(n+1,length(CT); X=a:h:b; Y(1,:)= CT;for k=1:nk1=feval(dydx,X(k),Y(k,:) x2=X(k)+h/2;y2=Y(k,:)+k1*h/2;k2=feval(dydx,x2,y2); k3=feval(dydx,x2,Y(k,:)+k2*h/2); k4=feval(dydx, X(k)+h,Y(k,:)+k3*h); Y(k+1,:)=Y(k,:)+h*(k1+2*k2+2*k3+k4)/6;k=k+1;endfor k=2:n+1wucha(k)=norm(Y(k)-Y(k-1); k=k+1;endX=X(1:n+1);Y=Y(1:n+1,:);k=1:n+1;wucha=wucha(1:k,:);P=k,X,Y,wucha;调用dydx.m求解，在MATLAB工作窗口输入程序CT=1;1;0;h=0.0001;k,X,Y,wucha,P=RK4z(dydx,0,0.01,CT,h),H=0.0001,0.01; x,y=ode15s(dydx,H,CT);plot(X,Y(:,1),g-,x,y(:,1),bo,X,Y(:,2),m:,x,y(:,2),cp,X,Y(:,3),r.,x,y(:,3),kd)xlabel(轴it x); ylabel(轴it y)title(分别用自定义函数和ode15s函数求解刚性方程方程组的图形)legend(用RK4z函数解刚性方程的y1的曲线,用ode15s函数解刚性方程的y1的曲线,用RK4z函数解刚性方程的y2的曲线,用ode15s函数解刚性方程的y2的曲线, 用RK4z函数解刚性方程的y3的曲线,用ode15s函数解刚性方程的y3的曲线)5.3用改进Euler方法的Matlab编程解法function x,y=eulerpro(dydx,a,b,CT,h);n=fix(b-a)/h);x=zeros(n+1,1); y=zeros(n+1,length(CT); y(1,:)= CT;for i=1:nx(i+1)=x(i)+h;f1=feval(dydx,x(i),y(i,:);y1(i+1,:)=y(i,:)+h*f1;f2=feval(dydx,x(i+1),y1(i+1,:);y2(i+1,:)=y(i,:)+h*f2;y(i+1,:)=(y1(i+1,:)+y2(i+1,:)/2;end调用dydx.m求解，在MATLAB工作窗口输入程序CT=1;1;0;h=0.0001;a=0;b=0.01;X,Y=eulerpro(dydx,a,b,CT,h),xlabel(轴it x); ylabel(轴it y)H=0.0001,0.01; x,y=ode15s(dy123,H,CT);plot(X,Y(:,1),g-,x,y(:,1),bo,X,Y(:,2),m:,x,y(:,2),cp,X,Y(:,3),r.,x,y(:,3),kd)xlabel(轴it x); ylabel(轴it y)title(分别用自定义函数和ode15s函数求解刚性方程方程组的图形)legend(用eulerpro函数解刚性方程的y1的曲线,用ode15s函数解刚性方程的y1的曲线,用eulerpro函数解刚性方程的y2的曲线,用ode15s函数解刚性方程的y2的曲线, 用eulerpro函数解刚性方程的y3的曲线,用ode15s函数解刚性方程的y3的曲线)6计算结果运行结果详见附录6.1用四阶Runge-Kutta方法的Matlab编程解法的结果以及与精确解的比较6.2用改进Euler方法的Matlab编程解法的结果以及与精确解的比较（注：计算结果详见附录1，附录2）7.结果分析由图可知次方法与精确解的契合度非常好，基本上与精确解保持一致，由此可见四阶Runge-Kutta方法是一种高精度的单步方法。此法对比改进Euler方法精确度更高。但是，相对的计算步骤比Euler改进方法要繁琐。综上，当计算低精度问题时可以使用改进Euler方法来处理问题，而如果精度要求较高，就要使用四阶Runge-Kutta方法。致谢本次课程设计主要针对一阶常微分方程组的初值问题，利用改进Euler方法以及Runge-Kutta方法解决一阶产微分方程的算法进行推广到一阶常微分方程组的算法。其中一部分编程思路参考了相关的参考书，经过反复验证然后得到的结果。其中在编写改进Euler方法的程序时，由于参考程序是非方程组的程序，在改写的时候忽略了向量的方向导致不匹配，通过仔细检查后成功运行程序。回顾起此次课程设计，至今我仍感慨颇多，的确，从找参考，设计到定稿，从理论到实践，在整整一星期的日子里，可以说得是苦多于甜，但是可以学到很多很多的的东西，同时不仅可以巩固以前所学过的知识，而且学到了很多在书本上所没有学到过的知识。通过这次课程设计使我懂得了理论与实际相结合是很重要的，只有理论知识是远远不够的，只有把所学的理论知识与实践相结合起来，从理论中得出结论，才能提高自己的实际动手能力和独立思考的能力。同时在设计的过程中发现了自己的不足之处，对以前所学过的知识理解得不够深刻，掌握得不够牢固。通过这次课程设计之后，一定把以前所学过的知识仔细复习，特别是要注意细节部分。参考文献1 李庆扬,王能超，易大义.数值分析M.清华大学出版社,2008. 2 宋叶志，贾东永.Matlab数值分析与应用M.机械工业出版社,2009. 3 黄明游,冯果忱.数值分析M.高等教育出版社,2008.附录改进Euler方法计算结果Y = 1.0000 1.0000 0 0.9050 1.0125 -0.0825 0.8172 1.0427 -0.1401 0.7351 1.0865 -0.1785 0.6579 1.1405 -0.2016 0.5855 1.2020 -0.2125 0.5177 1.2686 -0.2137 0.4547 1.3378 -0.2075 0.3965 1.4075 -0.1959 0.3435 1.4757 -0.1808 0.2956 1.5408 -0.1636 0.2527 1.6016 -0.1456 0.2149 1.6573 -0.1278 0.1817 1.7074 -0.1108 0.1530 1.7518 -0.0952 0.1283 1.7906 -0.0811 0.1072 1.8243 -0.0686 0.0893 1.8531 -0.0576 0.0742 1.8776 -0.0482 0.0615 1.8983 -0.0402 0.0509 1.9157 -0.0334 0.0420 1.9302 -0.0277 0.0347 1.9424 -0.0229 0.0286 1.9525 -0.0189 0.0235 1.9608 -0.0156 0.0194 1.9678 -0.0129 0.0159 1.9735 -0.0106 0.0131 1.9782 -0.0087 0.0108 1.9821 -0.0072 0.0088 1.9853 -0.0059 0.0073 1.9879 -0.0048 0.0060 1.9901 -0.0040 0.0049 1.9918 -0.0033 0.0040 1.9933 -0.0027 0.0033 1.9945 -0.0022 0.0027 1.9955 -0.0018 0.0022 1.9963 -0.0015 0.0018 1.9970 -0.0012 0.0015 1.9975 -0.0010 0.0012 1.9979 -0.0008 0.0010 1.9983 -0.0007 0.0008 1.9986 -0.0005 0.0007 1.9989 -0.0004 0.0006 1.9991 -0.0004 0.0005 1.9992 -0.0003 0.0004 1.9994 -0.0002 0.0003 1.9995 -0.0002 0.0002 1.9996 -0.0002 0.0002 1.9996 -0.0001 0.0002 1.9997 -0.0001 0.0001 1.9998 -0.0001 0.0001 1.9998 -0.0001 0.0001 1.9998 -0.0001 0.0001 1.9999 -0.0001 0.0001 1.9999 -0.0000 0.0001 1.9999 -0.0000 0.0000 1.9999 -0.0000 0.0000 1.9999 -0.0000 0.0000 1.9999 -0.0000 0.0000 1.9999 -0.0000 0.0000 1.9999 -0.0000 0.0000 2.0000 -0.0000 0.0000 2.0000 -0.0000四阶Runge-Kutta方法结算结果P = 1.0000 0 1.0000 1.0000 0 0 2.0000 0.0001 0.9045 1.0112 -0.0843 0.0955 3.0000 0.0002 0.8164 1.0401 -0.1435 0.0881 4.0000 0.0003 0.7343 1.0814 -0.1843 0.0821 5.0000 0.0004 0.6575 1.1311 -0.2114 0.0769 6.0000 0.0005 0.5855 1.1863 -0.2282 0.0719 7.0000 0.0006 0.5185 1.2445 -0.2369 0.0670 8.0000 0.0007 0.4565 1.3042 -0.2393 0.0620 9.0000 0.0008 0.3995 1.3638 -0.2367 0.0570 10.0000 0.0009 0.3475 1.4222 -0.2302 0.0520 11.0000 0.0010 0.3006 1.4786 -0.2207 0.0469 12.0000 0.0011 0.2586 1.5324 -0.2090 0.0420 13.0000 0.0012 0.2213 1.5830 -0.1957 0.0373 14.0000 0.0013 0.1884 1.6302 -0.1814 0.0328 15.0000 0.0014 0.1597 1.6737 -0.1666 0.0287 16.0000 0.0015 0.1348 1.7134 -0.1517 0.0249 17.0000 0.0016 0.1134 1.7495 -0.1370 0.0214 18.0000 0.0017 0.0950 1.7820 -0.1229 0.0184 19.0000 0.0018 0.0794 1.8111 -0.1095 0.0156 20.0000 0.0019 0.0662 1.8368 -0.0970 0.0132 21.0000 0.0020 0.0550 1.8596 -0.0854 0.0112 22.0000 0.0021 0.0456 1.8796 -0.0747 0.0094 23.0000 0.0022 0.0378 1.8971 -0.0651 0.0079 24.0000 0.0023 0.0312 1.9123 -0.0565 0.0066 25.0000 0.0024 0.0258 1.9254 -0.0488 0.0055 26.0000 0.0025 0.0212 1.9367 -0.0420 0.0045 27.0000 0.0026 0.0175 1.9465 -0.0360 0.0037 28.0000 0.0027 0.0144 1.9548 -0.0308 0.0031 29.0000 0.0028 0.0118 1.9619 -0.0262 0.0026 30.0000 0.0029 0.0097 1.9680 -0.0223 0.0021 31.0000 0.0030 0.0080 1.9731 -0.0189 0.0017 32.0000 0.0031 0.0066 1.9775 -0.0160 0.0014 33.0000 0.0032 0.0054 1.9811 -0.0135 0.0012 34.0000 0.0033 0.0044 1.9842 -0.0113 0.0010 35.0000 0.0034 0.0036 1.9868 -0.0095 0.0008 36.0000 0.0035 0.0030 1.9890 -0.0080 0.0007 37.0000 0.0036 0.0024 1.9908 -0.0067 0.0005 38.0000 0.0037 0.0020 1.9924 -0.0056 0.0004 39.0000 0.0038 0.0016 1.9937 -0.0047 0.0004 40.0000 0.0039 0.0013 1.9947 -0.0039 0.0003 41.0000 0.0040 0.0011 1.9956 -0.0033 0.0002 42.0000 0.0041 0.0009 1.9964 -0.0027 0.0002 43.0000 0.0042 0.0007 1.9970 -0.0023 0.0002 44.0000 0.0043 0.0006 1.9975 -0.0019 0.0001 45.0000 0.0044 0.0005 1.9979 -0.0016 0.0001 46.0000 0.0045 0.0004 1.9983 -0.0013 0.0001 47.0000 0.0046 0.0003 1.9986 -0.0011 0.0001 48.0000 0.0047 0.0003 1.9988 -0.0009 0.0001 49.0000 0.0048 0.0002 1.9990 -0.0007 0.0000 50.0000 0.0049 0.0002 1.9992 -0.0006 0.0000 51.0000 0.0050 0.0001 1.9993 -0.0005 0.0000 52.0000 0.0051 0.0001 1.9994 -0.0004 0.0000 53.0000 0.0052 0.0001 1.9995 -0.0003 0.0000 54.0000 0.0053 0.0001 1.9996 -0.0003 0.0000 55.0000 0.0054 0.0001 1.9997 -0.0002 0.0000 56.0000 0.0055 0.0001 1.9997 -0.0002 0.0000 57.0000 0.0056 0.0000 1.9998 -0.0002 0.0000 58.0000 0.0057 0.0000 1.9998 -0.0001 0.0000 59.0000 0.0058 0.0000 1.9998 -0.0001 0.0000 60.0000 0.0059 0.0000 1.9999 -0.0001 0.0000 61.0000 0.0060 0.0000 1.9999 -0.0001 0.0000 62.0000 0.0061 0.0000 1.9999 -0.0001 0.0000 63.0000 0.0062 0.0000 1.