我正在做以下操作来生成具有特定范围内的特征值的随机矩阵:
function mat = randEig(dim, rReal)
D=diff(rReal).*rand(dim,1)+rReal(1);
P=rand(dim);
mat=P*diag(D)/P;
end
但我也希望能够生成具有复杂(共轭)特征值的随机实矩阵。如何做到这一点?相似变换技巧将返回复杂矩阵。
编辑:好的,我设法通过搭载MATLAB的cdf2rdf
函数(基本上是下面的第二个函数)来做到这一点。
function mat = randEig(dim, rangeEig, nComplex)
if 2*nComplex > dim
error('Cannot happen');
end
if nComplex
cMat=diff(rangeEig).*rand(dim-2*nComplex,1)+rangeEig(1);
for k=1:nComplex
rpart=(diff(rangeEig).*rand(1,1)+rangeEig(1))*ones(2,1);
ipart=(diff(rangeEig).*rand(1,1)+rangeEig(1))*i;
ipart=[ipart; -ipart];
cMat=[cMat; rpart+ipart];
end
else
cMat=diff(rangeEig).*rand(dim,1)+rangeEig(1);
end
D=cMat;
realDform = comp2rdf(diag(D));
P=rand(dim);
mat=P*realDform/P;
end
function dd = comp2rdf(d)
i = find(imag(diag(d))');
index = i(1:2:length(i));
if isempty(index)
dd=d;
else
if (max(index)==size(d,1)) | any(conj(d(index,index))~=d(index+1,index+1))
error(message('Complex conjugacy not satisfied'));
end
j = sqrt(-1);
t = eye(length(d));
twobytwo = [1 1;j -j];
for i=index
t(i:i+1,i:i+1) = twobytwo;
end
dd=t*d/t;
end
end
但代码很丑陋,主要是多次调用rand
的方式很烦人)。如果有人想发布一次调用rand
的答案,并设法做到这一点,我一定会接受并赞成。
我做了一个电话或两个电话:
function mat = randEig(dim, rangeEig, nComplex)
if 2*nComplex > dim
error('Cannot happen');
end
if nComplex
cMat=diff(rangeEig).*rand(2*nComplex,1)+rangeEig(1);
cPart=cMat(1:nComplex)*i;
cMat(1:nComplex)=[];
cPart=upsample(cPart,2);
cPart=cPart+circshift(-cPart,1);
cMat=upsample(cMat,2);
cMat=cMat+circshift(cMat,1);
cMat=cMat+cPart;
cMat=[diff(rangeEig).*rand(dim-2*nComplex,1)+rangeEig(1); cMat];
else
cMat=diff(rangeEig).*rand(dim,1)+rangeEig(1);
end
D=cMat;
realDform = comp2rdf(diag(D));
P=rand(dim);
mat=P*realDform/P;
end
function dd = comp2rdf(d)
i = find(imag(diag(d))');
index = i(1:2:length(i));
if isempty(index)
dd=d;
else
if (max(index)==size(d,1)) | any(conj(d(index,index))~=d(index+1,index+1))
error(message('Complex conjugacy not satisfied'));
end
j = sqrt(-1);
t = eye(length(d));
twobytwo = [1 1;j -j];
for i=index
t(i:i+1,i:i+1) = twobytwo;
end
dd=t*d/t;
end
end
如果有人可以拨打一个电话或更短/更优雅的代码,欢迎他们发布答案。