chenjob21ding.m

来自「用MATLAB实现用C代码可控性和客观性的判断和转移矩阵的求解等各种功能。」· M 代码 · 共 62 行

A=input('please enter the matrix A= ');
B=input('please enter the matrix B= ');
T=input('please enter the sampling time T= ');
E=input('please enter the set error E= ');
L=1;
a=0;
[m,n]=size(A);
[p,q]=size(B);
z1=zeros(n);
z=zeros(n);
for i=1:1:n
    for j=1:1:n
a=a+abs(A(i,j)); %B为A的范数
end
end
for i=0:1:1
    z=z+(A*T)^i/prod(1:i);
    z1=z1+(A*T)^i*T/prod(1:(i+1));
end
m1=det(z1);
m0=det(z);
while(1)
M=zeros(m,n);
 for i=0:1:L
 M=M+(A*T)^i/prod(1:i);
end
if (m0>det(M))   %求最小值
    m0=det(M);
end
e=a*T/(L+2);%计算e来确定是否符合精度,当小于1时前提成立
F=(a*T)^(L+1)/prod(1:(L+1))*(1/(1-e));
if ((e<1)&&(F<=E*m))
L
disp('The approximation of state transition matrix is')
M
break
else 
    L=L+1;
end
end
L=1;
while(1)
M1=zeros(m,n);
 for i=0:1:L
 M1=M1+((A*T)^i*T/prod(1:i)*(i+1));
end
if (m1>det(M1))   %求最小值
    m1=det(M1);
end
e=a*T/(L+2);%计算e来确定是否符合精度,当小于1时前提成立
F=(a*T)^(L+1)/prod(1:(L+1))*(1/(1-e));
if ((e<1)&&(F<=E*m1))
L
M1=M1*B;
M1
break
else 
    L=L+1;
end
end