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

《微分方程组的龙格库塔公式求解matlab版.pdf》由会员上传分享，免费在线阅读，更多相关内容在教育资源-天天文库。

1、微分方程组的龙格-库塔公式求解matlab版南京大学王寻1.一阶常微分方程组考虑方程组y'fx,y,z,yx0y0z'gx,y,z,zxz00其经典四阶龙格-库塔格式如下：对于n=0,1,2,...,计算hyyK2K2KKn1n12346hzzL2L2LLn1n12346其中K1fxn,yn,zn,L1gxn,yn,znhhK1hL1hhK1hL1K2fxn,yn,zn,L2gxn,yn,zn222

2、222Kfxh,yhK2,zhL2,Lgxh,yhK2,zhL23nnn3nnn222222Kfxh,yhK,zhL,Lgxh,yhK,zhL4nn3n34nn3n3下面给出经典四阶龙格-库塔格式解常微分方程组的matlab通用程序：%marunge4s.m%用途：4阶经典龙格库塔格式解常微分方程组y'=f(x,y),y(x0)=y0%格式：[x,y]=marunge4s(dyfun,xspan,y0,h)%dyfun为向量函数f(x,y),x

3、span为求解区间[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(

4、:,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求刚性微分方程组y'0.01y99.99z,y02,z'100z,z010.01x100x100x的数值解，其解析解为：yee,ze解：首先编写M函数dyfun.m%dyfun.mfunctionf=dyfun(t,y)f(1)=-0.01

5、*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

6、(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('解析解');可以看出数值解和解析解的运算结果误差很小：解析解21.510.5y(x)z(x)0-0.5-50050100150200250300350400450500数值解21.510.5y(x)z(x)0-0.5-50050100150200250300350400450500例2：考虑下面的Lorenz方程组dx

7、xydtdyxyxzdtdzxyzdt参数,,适当的取值会使得系统趋于混沌状态。取=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(:);再编写运行

8、函数：clearall;closeall;clc;[t,y]=marunge4s(@mafun,[0500],[012],0.005);plot(y(1,:),y(3,:),'r');得到如下图所示的结果：6050403020100-25-20-15-10-505101520252.高阶微分方程组对于高