矩阵块对角

Python

未完，sympy的通解不知怎么使用。请问是否有人了解这方面工作，谢谢。

import random
import numpy as np
import sympy

#A_1,A_2,A_3
A_1=np.array([[2,1,0,0],[1,2,0,0],[0,0,1,2],[0,0,2,1]])
A_2=np.array([[0,0,1,0],[0,0,0,1],[1,0,0,0],[0,1,0,0]])
A_3=np.array([[0,0,0,1],[0,0,1,0],[0,1,0,0],[1,0,0,0]])
E=np.array([[1,0,0,0],[0,1,0,0],[0,0,1,0],[0,0,0,1]])


#求解AX-XA=0
#等价于求解 （E@A-At@E)vec X=0
#求张量积E@A-At@E(对称的不用转置)
A_1_kron=np.kron(E,A_1)-np.kron(A_1,E)
A_2_kron=np.kron(E,A_2)-np.kron(A_2,E)
A_3_kron=np.kron(E,A_3)-np.kron(A_3,E)
#张量积矩阵化
A_matrix_1_kron = sympy.Matrix(A_1_kron)
A_matrix_2_kron = sympy.Matrix(A_2_kron)
A_matrix_3_kron = sympy.Matrix(A_3_kron)
#在列上合并张量积 变成求解AX=0的形式
A=np.vstack((A_1_kron,A_2_kron,A_3_kron))






#求解X
A_matrix = sympy.Matrix(A)
#X=sympy.symarray('x',16)


x_11,x_12,x_13,x_14=sympy.symbols('x_11,x_12,x_13,x_14')
x_21,x_22,x_23,x_24=sympy.symbols('x_21,x_22,x_23,x_24')
x_31,x_32,x_33,x_34=sympy.symbols('x_31,x_32,x_33,x_34')
x_41,x_42,x_43,x_44=sympy.symbols('x_41,x_42,x_43,x_44')

x_11=0.2057   #给x附初值
x_12=0.1525
x_13=0.3036
'''
x_11=random.random()  #给x附初值
x_12=random.random() 
x_13=random.random() 
'''

X=sympy.Matrix([x_11,x_12,x_13,x_14,x_21,x_22,x_23,x_24,x_31,x_32,x_33,x_34,x_41,x_42,x_43,x_44])
print c,y,A,X
X=sympy.solve(A_matrix*X)
print X
#结果[{x_14: 0.303600000000000, x_21: 0.152500000000000, x_32: 0.303600000000000, x_22: 0.205700000000000, x_42: 0.303600000000000, x_31: 0.303600000000000, x_34: 0.152500000000000, x_23: 0.303600000000000, x_43: 0.152500000000000, x_33: 0.205700000000000, x_41: 0.303600000000000, x_44: 0.205700000000000, x_24: 0.303600000000000}]
Z=[]
for key in X:
    print key,X[key]
    Z.append(X[key])
print Z


#求X特征向量
#X矩阵化
X=np.array([[x_11,x_12,x_13,x_14],[x_21,x_22,x_23,x_24],[x_31,x_32,x_33,x_34],[x_41,x_42,x_43,x_44]])
print X


'''
a,P=np.linalg.eig(X) ##特征值赋值给a，对应特征向量赋值给b     
print P                       ##特征向量

#求P转置
P_t=np.transpose(P) 

#最终结果
print "P*(A1)P=",P_t*A_1*P
print "P*A(2)P=",P_t*A_2*P
print "P*(A3)P=",P_t*A_3*P
'''

MATLAB

%生成A矩阵和对角阵
A_1=[2 1 0 0;1 2 0 0;0 0 1 2;0 0 2 1];
A_2=[0 0 1 0;0 0 0 1;1 0 0 0;0 1 0 0];
A_3=[0 0 0 1;0 0 1 0;0 1 0 0;1 0 0 0];
%得到矩阵A列数c_A
[r_A,c_A] = size(A_1) ;
E=eye(c_A);

%求解AX-XA=0
%等价于求解 （E@A-A'@E)vecX=0,其中@指kron积
%求张量积E@A-At@E
A_1_kron=kron(E,A_1)-kron((A_1)',E);
A_2_kron=kron(E,A_2)-kron((A_2)',E);
A_3_kron=kron(E,A_3)-kron((A_3)',E);
%按列合并
A_kron=[A_1_kron;A_2_kron;A_3_kron];
[r_A_kron,c_A_kron] = size(A_kron) ;
%求解X通解
r=rank(A_kron);
X=null(A_kron,r);
%求X特解，参数为随机数
for i=1:(c_A_kron-r)
X_e=X_e+rand*X(:,i);
end
%得到矩阵X_e列数c
[r_X,c_X] = size(X_e) ;
%将vecX整理成n*n矩阵
X=(X_e(1:c_A,1))'; %初始化X
for j=2:c_A
X=[X;(X_e(((j-1)*c_A+1):(j*c_A),1))'];
end
%求X的特征矩阵P
[P,D]=eig(X);

%A_1,A_2,A_3块对角化
dia_A_1=P'*A_1*P;
dia_A_2=P'*A_2*P;
dia_A_3=P'*A_3*P;

%输出结果
X
P
dia_A_1
dia_A_2
dia_A_3

MATLAB求通解是毫秒级的速度，Python都是读秒计，代码还复杂了不少。果然术业有专攻，下次不装逼用Python做数学问题了。