loadingX0=[0;0;0;0];%初始状态。
disp('矩阵A为:')
A=[-0.0507 -3.861 0 -32.2;-0.00117 -0.5164 1 0;-0.000129 1.4168 -0.4932 0;0 0 1 0] %求解2!3!4！5！. n2=factorial(2); n3=factorial(3); n4=factorial(4); n5=factorial(5); I=eye(4);
disp('取最高次为5时M的值为：')
M=I+A*t+(1/n2)*A^2*t^2+(1/n3)*A^3*t^3+(1/n4)*A^4*t^4+(1/n5)*A^5*t^5 %求解本系统的状态转移矩阵。 B=[0;-0.0717;-1.645;0]; XT=M*(X0+B) %求解系统的状态响应。 %画图
subplot(2,2,1)
ezplot(XT(1,1),[0,2])%画出x(t)d第一个分量X1(t),并把它显示在左上角。 subplot(2,2,2)
ezplot(XT(2,1),[0,2]) %画出x(t)d第二个分量X2(t),并把它显示在右上角。 subplot(2,2,3)
ezplot(XT(3,1),[0,2]) %画出x(t)d第三个分量X3(t),并把它显示在左下角。 subplot(2,2,4)
ezplot(XT(4,1),[0,2]) %画出x(t)d第四个分量X1(t),并把它显示在右下角。
运行结果为：
图1为系统的单位脉冲响应x（t）的四个状态的仿真波形图。
10080604020000.20.40.60.811.21.41.61.82tx1(t)的仿真波形-1.50-1-2-2-3-4-2.500.20.40.60.811.2x3(t)的仿真波形1.41.61.82t-500.20.40.60.811.21.41.61.82tx4(t)的仿真波形-3-400.20.40.60.811.21.4x2(t)的仿真波形1.61.82t0-1-2
图1 定义法求单位脉冲响应仿真波形
2、利用拉氏反变换法求eAt及仿真
eAt的求解及求系统单位脉冲响应的Matlab程序如下：
syms M s d1 t XT X0;
A=[-0.0507 -3.861 0 -32.2;-0.00117 -0.5164 1 0;-0.000129 1.4168 -0.4932 0;0 0 1 0];
disp('矩阵A的行列式如下:'); d1=det(A); I=eye(4);
disp('[sI-A]^(-1)为:'); B=(s*I-A); C=inv(B); digits(4) C=vpa(C)
disp('状态转移阵为'); D=ilaplace(C); digits(4); M=vpa(D) X0=[0;0;0;0];
B=[0;-0.0717;-1.645;0]; XT=M*(X0+B) %求解系统的状态响应。 %画图
subplot(2,2,1) %画出x(t)d第一个分量X1(t),并把它显示在左上角。 ezplot(XT(1,1),[0,2])
subplot(2,2,2) %画出x(t)d第二个分量X2(t),并把它显示在右上角。 ezplot(XT(2,1),[0,2])
subplot(2,2,3) %画出x(t)d第三个分量X3(t),并把它显示在左下角。 ezplot(XT(3,1),[0,2])
subplot(2,2,4) %画出x(t)d第一个分量X4(t),并把它显示在右下角。 ezplot(XT(4,1),[0,2])
运行结果为：
图2为系统的单位脉冲响应x（t）的四个状态的仿真波形图。
12010080604020000.20.40.60.81x1((t)1.21.41.61.82t0-0.5-1-1.5-2-2.5-300.20.40.60.81x2(t)1.21.41.61.82t-1.5-2-2.5-3-3.500.20.40.60.81x3(t)1.21.41.61.82t0-1-2-3-400.20.40.60.81x4(t)1.21.41.61.82t
图2 拉氏反变换法求单位脉冲响应仿真波形
3、利用凯莱—哈密顿定理求eAt及仿真
eAt的求解及求系统单位脉冲响应的Matlab程序如下：
clear clc
syms s t a M a0 a1 a2 a3 b0 b1 b2 b3 X0 B XT I=eye(4);
A=[-0.0507 -3.861 0 -32.2;-0.00117 -0.5164 1 0;-0.000129 1.4168 -0.4932 0;0 0 1 0]; [V,D]=eig(A);
b1=D(1,1);b2=D(2,2);b3=D(3,3);b4=D(4,4);
c=[1 b1 (b1)^2 (b1)^3;1 b2 (b2)^2 (b2)^3;1 b3 (b3)^2 (b3)^3;1 b4 (b4)^2 (b4)^3]; d=inv(c);
disp('系数矩阵如下:');
f=[exp((b1)*t);exp((b2)*t);exp((b3)*t);exp((b4)*t)]; a=d*f;
a0=a(1,1); %求解公式中的第一个系数 a1=a(2,1); %求解公式中的第二个系数 a2=a(3,1); %求解公式中的第三个系数 a3=a(4,1); %求解公式中的第四个系数 disp('状态转移阵位');
M=a0*I+a1*A+a2*A^2+a3*A^3; %求系统的状态方程 X0=[0;0;0;0];
B=[0;-0.0717;-1.645;0]; XT=M*(X0+B)
subplot(2,2,1) %画出x(t)d第一个分量X1(t),并把它显示在左上角。 ezplot(XT(1,1),[0,2])
subplot(2,2,2) %画出x(t)d第二个分量X2(t),并把它显示在右上角。 ezplot(XT(2,1),[0,2])
subplot(2,2,3) %画出x(t)d第三个分量X3(t),并把它显示在左下角。 ezplot(XT(3,1),[0,2])
subplot(2,2,4) %画出x(t)d第一个分量X4(t),并把它显示在右下角。 ezplot(XT(4,1),[0,2])
运行结果为：
图3为系统的单位脉冲响应x（t）的四个状态的仿真波形图。
12010080604020000.20.40.60.81x1(t)1.21.41.61.82t0-0.5-1-1.5-2-2.5-300.20.40.60.81x2(t)1.21.41.61.82t-1.5-2-2.5-3-3.500.20.40.60.811.21.41.61.82t0-1-2-3-400.20.40.60.811.21.41.61.82tx3(t) x4(t)
图3 凯莱—哈密顿定理求单位脉冲响应仿真波形
4、利用约旦标准型法求eAt及仿真
eAt的求解及求系统单位脉冲响应的Matlab程序如下：
clear clc
syms M B N t a0 a1 a2 a3 X0 XT
A=[-0.0507 -3.861 0 -32.2;-0.00117 -0.5164 1 0;-0.000129 1.4168 -0.4932 0;0 0 1 0];
[T,D]=eig(A);%求特征值和特征向量. H=inv(T)*A*T; a0=D(1,1); a1=D(2,2); a2=D(3,3); a3=D(4,4);
N=[exp(a0*t) 0 0 0;0 exp(a1*t) 0 0;0 0 exp(a2*t) 0;0 0 0 exp(a3*t)]; M=T*N*inv(T); X0=[0;0;0;0];
B=[0;-0.0717;-1.645;0]; XT=M*(X0+B)
%其中，M即为所求的结果，B为A矩阵的约旦标准型，T为变换矩阵。 subplot(2,2,1) %画出x(t)d第一个分量X1(t),并把它显示在左上角。 ezplot(XT(1,1),[0,2])
subplot(2,2,2) %画出x(t)d第二个分量X2(t),并把它显示在右上角。 ezplot(XT(2,1),[0,2])
subplot(2,2,3) %画出x(t)d第三个分量X3(t),并把它显示在左下角。 ezplot(XT(3,1),[0,2])
subplot(2,2,4) %画出x(t)d第一个分量X4(t),并把它显示在右下角。 ezplot(XT(4,1),[0,2])
运行结果为：