微分方程组的龙格库塔公式求解matlab版

微分方程组的龙格-库塔公式求解matlab版南京大学王寻1.一阶常微分方程组考虑方程组0000zxz,z,y,xg'zyxy,z,y,xf'y其经典四阶龙格-库塔格式如下：对于n=0,1,2,...,计算4321143211226226LLLLhzzKKKKhyynnnn其中33433422322311211211222222222222hLz,hKy,hxgL,hLz,hKy,hxfKhLz,hKy,hxgL,hLz,hKy,hxfKhLz,hKy,hxgL,hLz,hKy,hxfKz,y,xgL,z,y,xfKnnnnnnnnnnnnnnnnnnnnnnnn下面给出经典四阶龙格-库塔格式解常微分方程组的matlab通用程序：%marunge4s.m%用途：4阶经典龙格库塔格式解常微分方程组y'=f(x,y),y(x0)=y0%格式：[x,y]=marunge4s(dyfun,xspan,y0,h)%dyfun为向量函数f(x,y),xspan为求解区间[x0,xn],%y0为初值向量，h为步长，x返回节点，y返回数值解向量function[x,y]=marunge4s(dyfun,xspan,y0,h)x=xspan(1):h:xspan(2);y=zeros(length(y0),length(x));y(:,1)=y0(:);forn=1:(length(x)-1)k1=feval(dyfun,x(n),y(:,n));k2=feval(dyfun,x(n)+h/2,y(:,n)+h/2*k1);k3=feval(dyfun,x(n)+h/2,y(:,n)+h/2*k2);k4=feval(dyfun,x(n+1),y(:,n)+h*k3);y(:,n+1)=y(:,n)+(h/6).*(k1+2*k2+3*k3+k4);end如下为例题：例1：取h=0.02，利用程序marunge4s.m求刚性微分方程组10100209999010z,z'z,y,z.y.'y的数值解，其解析解为：xxx.ez,eey100100010解：首先编写M函数dyfun.m%dyfun.mfunctionf=dyfun(t,y)f(1)=-0.01*y(1)-99.99*y(2);f(2)=-100*y(2);f=f(:);然后编写一个执行函数：closeall;clearall;clc;[x,y]=marunge4s(@dyfun,[0500],[21],0.002);plot(x,y);axis([-50500-0.52]);text(120,0.4,'y(x)');text(70,0.1,'z(x)');title('数值解');figure;y1=exp(-0.01*x)+exp(-100*x);%解析解，用来做对比的。z1=exp(-100*x);plot(x,y1,'r');holdon;plot(x,z1,'g');text(120,0.4,'y(x)');text(70,0.1,'z(x)');axis([-50500-0.52]);title('解析解');可以看出数值解和解析解的运算结果误差很小：-50050100150200250300350400450500-0.500.511.52y(x)z(x)解析解-50050100150200250300350400450500-0.500.511.52y(x)z(x)数值解例2：考虑下面的Lorenz方程组zxydtdzxzyxdtdyyxdtdx参数,,适当的取值会使得系统趋于混沌状态。取=30,=2.8,=12，利用经典四阶龙格-库塔法求其数值解，并绘制z随x变化的曲线。解：首先编写函数的m文件：%mafun.mfunctionff=mafun(t,y)b=2.8;r=30;sigma=12;ff(1)=-sigma*y(1)+sigma*y(2);ff(2)=r*y(1)-y(2)-y(1)*y(3);ff(3)=y(1)*y(2)-b*y(3);ff=ff(:);再编写运行函数：clearall;closeall;clc;[t,y]=marunge4s(@mafun,[0500],[012],0.005);plot(y(1,:),y(3,:),'r');得到如下图所示的结果：-25-20-15-10-5051015202501020304050602.高阶微分方程组对于高阶微分方程，总是可以化成方程组的形式。例如，二阶方程0000'yx'y,yxy'y,y,xgy总是可以化为一阶方程组：00000z'yxz,yxyz,y,xg'zz'y因此不需要对于高阶方程给出计算公式。例3：求二阶方程11151123'yy.x,yy的数值解，其解析解为：21xy解：首先将二阶方程写为一阶方程组的形式：112113z,y'zy,z'y再编写函数的m文件：%这是一个二阶微分方程functionf=dyfun1(t,y)f(1)=y(2);%y(1)是y,y(2)是z。f(2)=2*(y(1)*y(1)*y(1));f=f(:);直接使用ode45来计算：closeall;clearall;clc;[x,y]=ode45('dyfun1',[11.5]',[-1-1]')plot(x,(y(:,1)));其结果如图：11.051.11.151.21.251.31.351.41.451.5-2.2-2-1.8-1.6-1.4-1.2-1-0.8这与解析解作出来的图像很接近。