![]() |
#2
jinanzht2009-11-17 22:10
|
程序:
function rk(A,x,h,y0)
i=1;
y(i,:)=y0;
m=length(A);
b=x(length(x));
while x(i)<b
a=[x(i),y(i,:)];
for l=1:m
k1(l)=eval(A{l},a);
end
a=[x(i)+h/2,y(i,:)+h/2*k1];
for l=1:l
k2(l)=eval(A{l},a);
end
a=[x(i)+h/2,y(i,:)+h/2*k2];
for l=1:m
k3(l)=eval(A{l},a);
end
a=[x(i)+h,y(i,:)+h*k3];
for l=1:m
k4(l)=eval(A{l},a);
end
y(i+1,:)=y(i,:)+h/6*(k1+2*k2+2*k3+k4);
i=i+1;
end
y,k1,k2,k3,k4
运行结果:
y1='-10*a(2)-9*a(3)+8';
y2='-9*a(2)-8*a(3)+7';
A=[{y1},{y2}];
rk(A,0:0.05:1,0.05,ones(1,2))
y =
1.0000 1.0000
0.6417 0.6716
0.4972 0.5345
0.4403 0.4758
0.4193 0.4491
0.4129 0.4356
0.4126 0.4275
0.4148 0.4216
0.4181 0.4166
0.4217 0.4120
0.4255 0.4075
0.4294 0.4030
0.4334 0.3986
0.4373 0.3941
0.4413 0.3897
0.4453 0.3852
0.4493 0.3807
0.4533 0.3762
0.4574 0.3717
0.4614 0.3672
0.4655 0.3627
k1 =
0.0809 -0.0904
k2 =
0.0811 -0.0906
k3 =
0.0811 -0.0906
k4 =
0.0812 -0.0907
我想把中间的K1,k2,k3,k4都显示出来,该怎么做啊,请教赐教