4阶龙格库塔法的matlab仿真程序

a=[010;001;-22.06-27-10];
b=[0;0;40.6];
u=1;
x=[0;0;0];
h=0.4455;
y=x‘;
t（1）=0;
for i=1:100
    t（i+1）=h*i;
    k1=a*x+b*u;
    k2=a*（x+h/2*k1）+b*u;
    k3=a*（x+h/2*k2）+b*u;
    k4=a*（x+h*k3）+b*u;
    x=x+h/6*（k1+2*k2+2*k3+k4）;
    y=[y;x‘];
    y3=y‘;
    y2=y3（1:）;
end

subplot（221）;
plot（ty2）;
hold on;
y1=1.84-4.95*t.*exp（-1.88*t）-1.5*exp（-1.88*t）-0.34*exp（-6.24*t）;
plot（ty1‘r‘）;
grid
xlabel（‘time（sec）‘）;
ylabel（‘y‘）;
legend（‘y1解析‘‘y2数值‘）;
title（‘龙格——库塔法 临界步长h=0.4455‘）;

a=[010;001;-22.06-27-10];
b=[0;0;40.6];
u=1;
x=[0;0;0];
h=0.22275;
y=x‘;
t（1）=0;
for i=1:100
    t（i+1）=h*i;
    k1=a*x+b*u;
    k2=a*（x+h/2*k1）+b*u;
    k3=a*（x+h/2*k2）+b*u;
    k4=a*（x+h*k3）+b*u;
    x=x+h/6*（k1+2*k2+2*k3+k4）;
    y=[y;x‘];
    y3=y‘;
    y2=y3（1:）;
end

subplot（222）;
plot（ty2）;
hold on;
y1=1.84-4.95*t.*exp（-1.88*t）-1.5*exp（-1.88*t）-0.34*exp（-6.24*t）;
plot（ty1‘r‘）;
grid
xlabel（‘time（sec）‘）;
ylabel（‘y‘）;
legend（‘y1解析‘‘y2数值‘）;
title（‘龙格——库塔法 1/2临界步长h=0.22275‘