龙格库塔方法matlab实现

龙格库塔方法matlab实现~

function ff=rk(yy,x0,y0,h,a,b)%yy为y的导函数,x0,y0,为初值,h为步长,a,b为区间

c=(b-a)/h+1;i1=1; %c为迭代步数；i1为迭代步数累加值

y=y0;z=zeros(c,6); %z生成c行，5列的零矩阵存放结果；

%每行存放c次迭代结果，每列分别存放k1~k4及y的结果

for x=a:h:b

if i1<=c

k1=feval(yy,x,y);

k2=feval(yy,x+h/2,y+(h*k1)/2);

k3=feval(yy,x+h/2,y+(h*k2)/2);

k4=feval(yy,x+h,y+h*k3);

y=y+(h/6)*(k1+2*k2+2*k3+k4);

z(i1,1)=x;z(i1,2)=k1;z(i1,3)=k2;z(i1,4)=k3;z(i1,5)=k4;z(i1,6)=y; i1=i1+1;

end

end

fprintf(‘结果矩阵，第一列为x(n),第二列~第五列为k1~k4，第六列为y（n+1）的结果') z

%在命令框输入下列语句

%yy=inline('x+y');

%>> rk(yy,0,1,0.2,0,1)

%将得到结果

%结果矩阵，第一列为x(n),第二列~第五列为k1~k4第六列为y（n+1）的结果

%z =

% 0 1.0000 1.2000 1.2200 1.4440 1.2428

% 0.2000 1.4428 1.6871 1.7115 1.9851 1.5836

% 0.4000 1.9836 2.2820 2.3118 2.6460 2.0442

% 0.6000 2.6442 3.0086 3.0451 3.4532 2.6510

% 0.8000 3.4510 3.8961 3.9407 4.4392 3.4365

% 1.0000 4.4365 4.9802 5.0345 5.6434 4.4401