%--------------------------------------------------------------------------
% process for subspace iteration method
% A0: The initial modal matrix
% eps: Accuracy value
% num: the number of iterations
% n: Number of DOFs
% s: Iterative order
% D: Dynamic matrix
% W: The square of Natural frequency
% A: Main modes
% by: ZhouWeiyan 2017/11/28
%--------------------------------------------------------------------------
eps = 1e-4;
num = 1000;
n = size(K,1);
itetimes = 0; %用于保存迭代的次数
s=p; % 迭代前s阶,快速收敛到所求前p阶
%% Input initial martix A0
A=[1 2 3 1 0;1 2 -1 0 1]';
Ac=[0.3621 0.6956 0.8450 0.9212 1;0.7726 1 0.2108 -0.2925 -0.9937]'; %accurate the first two modes
N2Err1=zeros(1,num); %N2Err1 is used to save the error of the first mode in the iteration process
N2Err2=zeros(1,num); %N2Err2 is used to save the error of the second mode in the iteration process
%% subspace iteration method
D = inv(K)*M; % D = M\K; dynamic matrix
for j = 1:num
Aj = D*A;
Aj = column_normalization(Aj); %matrix column elements normalized
Mj = Aj'*M*Aj;
Kj = Aj'*K*Aj;
[aj,W] = eig(Mj\Kj);
Ajj = column_normalization(Aj*aj); %matrix column elements normalized
N2Err1(j)=norm((Ajj(1,:)-Ac(1,:)));
N2Err2(j)=norm((Ajj(2,:)-Ac(2,:)));
if( norm(Ajj-Ac) <= eps) %meet precision requirement
disp(['Ending at ',num2str(j),' th iteration']);
A = Ajj;
itetimes=j;
break;
else
A = Ajj;
end
W = eig(Mj\Kj);
[W_sort W] = sort(W,'ascend');% Reorder.%other one is 'descend'
%排序,W_sort为排序后特征值向量
A = A(:,W);
end
%% output the resualt
%disp(W_sort);
for i = 1 : p
W_sort(i) = sqrt(W_sort(i));
AAi = A(:,i);
disp([num2str(i),' th modal found with frequency as ',num2str(W_sort(i)),' Hz',' ...']);
disp([num2str(i),' th modal found with mode as ']);
disp(AAi);
end
%% plot the result
N2Err1_used=zeros(1,itetimes);
N2Err2_used=zeros(1,itetimes);
for i=1:itetimes
N2Err1_used(i)=N2Err1(i);
N2Err2_used(i)=N2Err2(i);
end
x=1:itetimes;
plot(x,N2Err1_used);
hold on
plot(x,N2Err2_used,'r:');
grid on
ylabel('两阶特征值的各次迭代值与准确解的二范数之差');
xlabel('迭代次数');
legend('第一阶特征值','第二阶特征值',1);
评论5