matlab求齐次线性微分方程组的基解矩阵:某过程求解释

Question

A=[2,0,0,0,0;4,2,4,0,-4;4,-1,6,1,-5;12,-2,12,4,-14;8,-1,8,1,-7]
[v,lambda]=eig(sym(A));n=length(lambda);
if length(v(1,:))==n syms t
Phi=[];for j=1:n
Phi=[Phi,eval(v(:,j)*exp(lambda(j,j)*t))];
end;Phi=Phi
Psi=eval(Phi/subs(Phi,t,'0'))
else disp('False');end

请说一下上面程序（最好各个步骤）的意思~是关于求x'=Ax的基解矩阵的~
A =
2 0 0 0 0
4 2 4 0 -4
4 -1 6 1 -5
12 -2 12 4 -14
8 -1 8 1 -7

Answer 1

A=[2,0,0,0,0;4,2,4,0,-4;4,-1,6,1,-5;12,-2,12,4,-14;8,-1,8,1,-7]

% 下面这句求矩阵A的特征值和特征向量，以特征值构成对角阵
lambda[v,lambda]=eig(sym(A));
% n本质上是系统的阶次
n=length(lambda);

% 这个判断有点莫名其妙，我想不出什么情况下if的条件会不满足
if length(v(1,:))==n
syms t % 定义符号变量t
Phi=[]; % 置初始的Φ(t)矩阵为空，后面依次添加各列
for j=1:n
% 通过循环依次添加Φ(t)矩阵各列，函数eval纯属多余
Phi=[Phi,eval(v(:,j)*exp(lambda(j,j)*t))];
end
% 下面这一句又是莫名其妙，如果想显示变量Phi，直接写变量名即可，没必要给自身赋值
Phi=Phi
% 通过处以Φ(0)，对Φ(t)进行类似归一化的操作。eval同样是多余的
Psi=eval(Phi/subs(Phi,t,'0'))
else
% 这个分支根本似乎不可能执行到
disp('False');
end

%==========================================
% 其实求状态转移矩阵更简单的做法如下：
[P, D]=eig(sym(A));
Psi1 = P * expm(D*t) * P^-1

% 可以验证，这样求出的基解矩阵和上面相同
Psi-Psi1

追问

我试过运行，当A=[2,1;0,2]的时候，就返回False...if条件不会满足，就应该是特征值为重根的时候？

追答

嗯，特征根有重根的情况确实会出现False。我之前有点想当然了。

这种情况下，矩阵应该按照约旦标准形分解，即状态转移矩阵为：
[P, D]=jordan(sym(A));
Psi1 = P * expm(D*t) * P^-1

另外，纠正一点，前面贴的代码有点小错误：
lambda[v,lambda]=eig(sym(A));
其中前面的lambda应属于上一行注释的内容。因为贴到知道里面经常会乱，修改时没注意。

Answer 2

小白酱油飘过。大神无视之。