Taylorss 发表于 2014-4-9 16:22

关于李雅普诺夫微分方程精细积分的一个程序,算不出数,求大....

M=;
C=;
K=;
In=ones(6);
g=1;
=eig(K,M);
A=;
N=20;
S0=142.75;
t=0.02;
it=t/2^N;
tt=it;
Ta=A*it+(A*it)^2*(In+(A*it)/3+(A*it)^2/12)/2;
fi=expm(A*tt);
G=;
Q0=G*(S0*ones(3))*G';
P0=Q0*it+(A*Q0+Q0*A')*it^2/2+(A^2*Q0+2*A*Q0*A'+Q0*(A^2)')*it^3/6+(A^3*Q0+3*A^2*Q0*A'+3*A*Q0*(A^2)'+Q0*(A^3)')*it^4/24;
P1=Q0*it^2/2+(A*Q0+Q0*A')*it^3/6+(A^2*Q0+2*A*Q0*A'+Q0*(A^2)')*it^4/24;
for tt=0:t/2^N:t/2
    P0=P0+eps;
    P1=P1+eps;
    P1=P1+(In+Ta*2+Ta*Ta)*P1*(In+Ta*2+Ta*Ta)'+2*tt*P0;
    P0=P0+(In+Ta*2+Ta*Ta)*P0*(In+Ta*2+Ta*Ta)';
end
for t2=t/2:t/2^N:t
      P1=P1+fi*fi*P1*(fi*fi)'+2*t2*P0;
      P0=P0+fi*fi*P0*(fi*fi)';
end
    Pk1=P1;
    Pk0=P0;
    for t1=t:t:N*t
      fi1=expm(A*t1);
      Pk1=fi1*Pk1*fi1'+t1*Pk0+Pk1;
      Pk0=fi1*Pk0*fi1'+Pk0;
      disp(Pk1)
      disp(Pk0)
    end

页: [1]
查看完整版本: 关于李雅普诺夫微分方程精细积分的一个程序,算不出数,求大....