《常微分方程组初值问题数值解的实现和算法分析》由会员分享，可在线阅读，更多相关《常微分方程组初值问题数值解的实现和算法分析（19页珍藏版）》请在金锄头文库上搜索。
1、常微分方程组初值问题数值解的实现和算法分析摘 要本次课程设计主要内容是用改进Euler方法和四阶Runge-Kutta方法解决常微分方程组初值问题的数值解法，通过分析给定题目使用Matlab编写程序计算结果并绘图然后区别两种方法的使用范围。最后对计算结果进行分析，得到结论。关键词：改进Euler，Runge-Kutta，初值问题目 录1前言12题目叙述13解题思路13.1一阶常微分方程的初值问题13.2一阶常微分方程组的初值问题23.2.1用Runge-Kutta方法计算解决一阶微分方程组初值问题的基本思路23.2.2用改进Euler方法计算解决一阶微分方程组初值问题的基本思路44用matla
2、b语言编程解决相关问题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前言常微分方程是解决工程实例的常用的工具，建立微分方程只是解决问题的第一步，通常需要求出方程的解来
3、说明实际现象，并加以检验。如果能得到解析形式的解固然是便于分析和应用的，但是我们知道，虽然求解常微分方程有各种各样的解析方法，但解析方法只能用来求解一些特殊类型的方程，实际问题终归结出来的微分方程主要靠数值解法1。数值解法就是一个十分重要的手段，而Euler方法以及Runge-Kutta方法又是数值解法最基础最常用的方法。通过对两种方法计算结果的对比，从而找到一个最为合适的计算方法。2题目叙述对以下的常微分方程组初值问题(1) 用四阶RungeKutta方法求解，根据计算结果画出解的图形(2) 用改进Euler方法求解，根据计算结果画出解的图形。3解题思路3.1一阶常微分方程的初值问题一阶常微
4、分方程的初值问题的一般形式为： 我们知道，只要函数f (x, y)适当光滑譬如关于y 满足利普希茨(Lipschitz)条件，即如果存在实数，使得理论上就可以保证初值问题的解存在并且唯一。所谓数值解法， 就是求上述问题在一系列离散点的近似值。两个相邻节点的间距称为步长13.2一阶常微分方程组的初值问题解一阶微分方程组时类似于单个方程的数值解法，只要把和理解为向量，那么就有：一阶微分方程组初值问题的形式为： （1）式中引入向量符号： 则（1）可写为：(2)3.2.1用Runge-Kutta方法计算解决一阶微分方程组初值问题的基本思路(2)式形式上与常数微分方程初值问题是一样的，只要注意向量函数运
5、算及其表示，就可以用初值问题的求解格式得到常微分方程组初值问题（2）的求解格式，由初值问题的经典Runge-kutta公式可得一阶常微分方程组初值问题（2）的Runge-kutta公式：注意上式是向量形式，其对应的分量形式为：微分方程理论告诉我们，高阶微分方程可转化为一阶微分方程组来研究，因此可以用一阶微分方程组初值问题揭发来解高阶微分方程初值问题。高阶微分方程初值问题的形式为： (3)令则（2）化为了一阶微分方程组初值问题：Runge-kutta方法巧妙利用函数在一些点上的函数值的线性组合，获得了高阶的数值解法，它避开了要获得高阶方法须对求高阶导数的不便，是离散化方法中Tayl情况，其中在准
6、确性的工作量的综合效果看，经典的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
7、,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+
8、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 输入计
9、算题目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:n
10、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;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,w
11、ucha;调用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函数解
12、刚性方程的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(dy
13、dx,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),
14、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方法。致谢本次课程设计主要针对一阶常微分方
