数值分析课程设计常微分方程组初值问题数值解的实现和算法分析毕业论文Word文档下载推荐.docx

1、（2）3.2.1用Runge-Kutta方法计算解决一阶微分方程组初值问题的基本思路（2）式形式上与常数微分方程初值问题是一样的，只要注意向量函数运算与其表示，就可以用初值问题的求解格式得到常微分方程组初值问题（2）的求解格式，由初值问题的经典Runge-kutta公式可得一阶常微分方程组初值问题（2）的Runge-kutta公式：注意上式是向量形式，其对应的分量形式为：微分方程理论告诉我们，高阶微分方程可转化为一阶微分方程组来研究，因此可以用一阶微分方程组初值问题揭发来解高阶微分方程初值问题。高阶微分方程初值问题的形式为：（3）令则（2）化为了一阶微分方程组初值问题：Runge-kutta方

2、法巧妙利用函数在一些点上的函数值的线性组合，获得了高阶的数值解法，它避开了要获得高阶方法须对求高阶导数的不便，是离散化方法中Tayl情况，其中在准确性的工作量的综合效果看，经典的Runge-kutta方法是首选or展开法的一个应用。Runge-kutta方法主要用于定步长的。Runge-kutta方法也常用于对多步法提供初值。3.2.2用改进Euler方法计算解决一阶微分方程组初值问题的基本思路改进Euler方法需要用Euler方法求出一个预测值然后再用梯形公式校正一次得到，即所求结果的迭代格式。 （4）为了方便编程可将（4）式改变为如下格式4用matlab语言编程解决相关问题4.1四阶Run

3、ge-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:x0+h/2）;k3=feval（f,k2*

4、（h/2）+T（1:k4=feval（f,k3*h+T（1:x0+h）;x0=x0+h;,t）=T（1:r-s,t-1）+（k1+k2*2+k3*2+k4）*（h/6）;End24.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）=

5、（y1+y2）/2;5 编程解决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

6、;Y（1,:）= CT;for k=1:k1=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;for k=2:wucha（k）=norm（Y（k）-Y（k-1）; k=k+1;X=X（1:n+1）;Y=Y（1:n+1,:）;k=1:n+1;wucha=wucha（1:k,:P=k,X

7、,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（:bo,X,Y（:,2）,m:cp,3）,r.kd）xlabel（轴it x ylabel（轴it ytitle（分别用自定义函数和ode15s函数求解刚性方程方程组的图形legend（用RK4z函数解刚性方程的y1的曲线,用ode15s函数解刚性方程的y1的曲线用RK4z函数解刚性方程的y2的曲线用od

8、e15s函数解刚性方程的y2的曲线, 用RK4z函数解刚性方程的y3的曲线用ode15s函数解刚性方程的y3的曲线5.3用改进Euler方法的Matlab编程解法function x,y=eulerpro（dydx,a,b,CT,h）;x=zeros（n+1,1）;y=zeros（n+1,length（CT）;y（1,:f1=feval（dydx,x（i）,y（i,:）;y1（i+1,:）=y（i,:）+h*f1f2=feval（dydx,x（i+1）,y1（i+1,:y2（i+1,:）+h*f2y（i+1,:）=（y1（i+1,:）+y2（i+1,:）/2;a=0;b=0.01;X,Y=eu

9、lerpro（dydx,a,b,CT,h）,dy123用eulerpro函数解刚性方程的y1的曲线用eulerpro函数解刚性方程的y2的曲线用eulerpro函数解刚性方程的y3的曲线6计算结果运行结果详见附录6.1用四阶Runge-Kutta方法的Matlab编程解法的结果以与与精确解的比较6.2用改进Euler方法的Matlab编程解法的结果以与与精确解的比较（注：计算结果详见附录1，附录2）7.结果分析由图可知次方法与精确解的契合度非常好，基本上与精确解保持一致，由此可见四阶Runge-Kutta方法是一种高精度的单步方法。此法对比改进Euler方法精确度更高。但是，相对的计算步骤比E

10、uler改进方法要繁琐。综上，当计算低精度问题时可以使用改进Euler方法来处理问题，而如果精度要求较高，就要使用四阶Runge-Kutta方法。致本次课程设计主要针对一阶常微分方程组的初值问题，利用改进Euler方法以与Runge-Kutta方法解决一阶产微分方程的算法进行推广到一阶常微分方程组的算法。其中一部分编程思路参考了相关的参考书，经过反复验证然后得到的结果。其中在编写改进Euler方法的程序时，由于参考程序是非方程组的程序，在改写的时候忽略了向量的方向导致不匹配，通过仔细检查后成功运行程序。回顾起此次课程设计，至今我仍感慨颇多，的确，从找参考，设计到定稿，从理论到实践，在整整一星期

11、的日子里，可以说得是苦多于甜，但是可以学到很多很多的的东西，同时不仅可以巩固以前所学过的知识，而且学到了很多在书本上所没有学到过的知识。通过这次课程设计使我懂得了理论与实际相结合是很重要的，只有理论知识是远远不够的，只有把所学的理论知识与实践相结合起来，从理论中得出结论，才能提高自己的实际动手能力和独立思考的能力。同时在设计的过程中发现了自己的不足之处，对以前所学过的知识理解得不够深刻，掌握得不够牢固。通过这次课程设计之后，一定把以前所学过的知识仔细复习，特别是要注意细节部分。参考文献1 庆扬,王能超，易大义.数值分析M.清华大学,2008. 2 宋叶志，贾东永.Matlab数值分析与应用M.

12、机械工业,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.6

13、016 -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.

14、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

15、.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.9999 -0.0001 0.0001 1

16、.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

17、 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.

18、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

19、.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

20、.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

21、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

22、 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

23、 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.000

24、3 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.999

25、9 -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.9999 -0.0000 0.0000 64.0000 0.0063 0.0000 1.9999 -0.0000 0.0000 65.0000 0.0064 0.0000 1.9999 -0.0000 0.0000 66.0000 0.0065 0.0000 1.9999 -0.0000 0.0000 67.0000 0.0066 0.0000 2.0000 -0.0000 0.0000 101.0000 0.0100 0.0000 2.0000 -0.0000 0.0000翻译关于feval的Matlab help原文如下： FEVAL Execute the specified function. FEVAL（F,x1,.,xn） evaluates the function specified by a function handle or function name,