%建立质量矩阵M，刚度矩阵K\
syms k m;
M=[m 0 0
0 m 0
0 0 2*m];
K=[2*k -k 0
-k 3*k -2*k
0 -2*k 2*k];
%*********迭代第n阶主阵型********************
%n为计数器
%
n=1;
while n<=3
%计算系统动力矩阵A
if n==1
A(:,:,n)=K\M;
elseif n~=1
A(:,:,n)=A(:,:,n-1)-(MP(:,n-1)\t(:,n-1))*f(:,n-1)*f(:,n-1)'*M;
end
%定义初始迭代向量X(1)
if n==1
X(:,1,n)=[1 1 1]';
elseif n==2
X(:,1,n)=[1 1 -1]';
elseif n==3
X(:,1,n)=[1 -1 1]';
end
%迭代过程，Y为中间矩阵,i为迭代次数
i=1;
Y(:,1,n)=A(:,:,n)*X(:,1,n);
X(:,2,n)=Y(3,1,n).\Y(:,1,n);
while abs(X(1,i,n)-X(1,i+1,n))>=0.000001&&abs(X(2,i,n)-X(2,i+1,n))>0.000001
Y(:,i+1,n)=A(:,:,n)*X(:,i+1,n);
X(:,i+2,n)=Y(3,i+1,n).\Y(:,i+1,n);
X(:,2);
i=i+1;
end
X(:,:,n)
f(:,n)=X(:,i,n);
t(n)=Y(3,i,n);
MP(:,n)=f(:,n)'*M*f(:,n);
n=n+1;
end
%输出数据过程
disp('第一阶主阵型的迭代结果');
X(:,:,1)
disp('第二阶主阵型的迭代结果');
X(:,:,2)
disp('第三阶主阵型的迭代结果');
X(:,:,3)
disp('φi的计算结果，矩阵的每列分别是1,2 ，3阶的');
f
disp('######******注意***********主阵型的迭代结果后面的0是系统的占位符号，不算计算结果');
