Matlabで時間応答（step, lsim, ode45, ode23s）

clear; close all;

g=32.174;%ヤード・ポンド法なので9.8ではない
Za=-236;
Zu=-0.227;
Xa=14.7;
Xu=-0.0215;
Madot=-0.28;
Ma=-3.78;
Mu=0;
Mq=-0.992;
Zde=-12.9;
Mde=-2.48;
U0=293.8;
Zq=-5.76;

global A; global B;
A=[Xu Xa 0 -g;
    Zu/U0,Za/U0 (U0+Zq)/U0 0;
    Mu+Madot*Zu/U0 Ma+Madot*Za/U0 Mq+Madot*(U0+Zq)/U0 0;
    0 0 1 0];   
B=[0; Zde/U0; Mde+Madot*Zde/U0; 0];
C=[1 0 0 0;
   0 1 0 0;
   0 0 1 0 ;
   0 0 0 1];
D=[0;0;0;0];
sys=ss(A,B,C,D);

opts = stepDataOptions('InputOffset',0,'StepAmplitude',0.1);

%step()を用いた時間応答
step(sys,opts);
xlim([0,50]);

%lsim()を用いた時間応答
N=200;
ts=linspace(0,50,N);
us=0.1*ones(N,1);
y=lsim(sys,us,ts);

figure();
subplot(411);
plot(ts,y(:,1));
subplot(412);
plot(ts,y(:,2));
subplot(413);
plot(ts,y(:,3));
subplot(414);
plot(ts,y(:,4));

%ode45を用いた時間応答
x0=[0 0 0 0 ]';
tspan=[0 50];
[t,x]=ode45(@dynamical_system,tspan,x0);
figure();
subplot(411);
plot(t,x(:,1));
subplot(412);
plot(t,x(:,2));
subplot(413);
plot(t,x(:,3));
subplot(414);
plot(t,x(:,4));

%ode23s()を用いた時間応答
[t,x]=ode23s(@dynamical_system,tspan,x0);
figure();
subplot(411);
plot(t,x(:,1));
subplot(412);
plot(t,x(:,2));
subplot(413);
plot(t,x(:,3));
subplot(414);
plot(t,x(:,4));

function dx = dynamical_system(t,x)
	% x = [u;α;q;θ; β;p;r;φ;ψ]
	global A; global B;
    delta_e=0.1;
    dx = A*x+B*delta_e;
end

４つの方法による計算は一致する