[振动力学] 课程考核报告：Matlab 实现邓克利法、瑞利法、里兹法、矩阵迭代法
作者：互联网
《振动理论》课程考核
1.建立系统柔度矩阵、刚度矩阵、质量矩阵，形成系统动力学方程；
2.求系统各阶真实固有频率和简正模态，并画出各阶简正模态的形状；
3.利用邓克利法求系统基频；
4.利用刚度矩阵表示的瑞利法求系统基频；
5.利用刚度矩阵表示的里兹法求系统前两阶固有频率；
6.利用动力矩阵表示的里兹法求系统前两阶固有频率；
7.利用矩阵迭代法求系统前两阶固有频率和简正模态；
8.比较 2、3、4（1）、5（2）和 6（2）计算基频的结果，列出五者的大小关系，作误差分析，说明现象；
9.总结 1 至 8 的计算实践，当中应说明各近似计算方法计算精度和收敛速度的影响因素；
代码：
% 作者：中山大学 2019 级 廉
% 第一问
% 符号定义
syms x1(t) x2(t) x3(t);
syms m E I L w lambda positive real;
% 向量的余弦近似度的计算公式
cosSim = @(a,b) sum(a.*b)/sqrt(sum(a.^2)*sum(b.^2));
% 无量纲化参数 dimenisonlessPara
dPara = (E*I/L^3)/m;
% 速度自由度 方便求质量矩阵 之后可能会被覆盖
x = [x1 x2 x3]';
dx = diff(x,1);
d2x = diff(x,2);
% 各点位置 方便求柔度矩阵 之后可能会被覆盖
x_pos = [L/4 L/2 3*L/4]';
% 求质量矩阵 M
M = m*[1 0 0;0 1 0;0 0 3];
% 求柔度矩阵 F
% 柔度矩阵是对称的，为了避免材料力学公式出现 xb > xa，可以先只求上三角
% 如果这里求错了，可能会使得求出来的 F 是奇异的
% 材料力学公式 长为 L 刚度 EI 的梁，B 点作用有单位力时 A 点的挠度
fab = @(a,b) a*b/(6*E*I*L)*(L^2-a^2-b^2);
% 柔度矩阵初始化
F = sym(zeros(size(x,1),size(x,1)));
% 由公式得柔度矩阵
for i = 1:1:3
for j = i:1:3
F(i,j) = fab(x_pos(i,1),x_pos(4-j,1));
end
end
% 复制上三角到下三角
for i = 2:1:3
for j = 1:1:i-1
F(i,j) = F(j,i);
end
end
% 求刚度矩阵 K
K = inv(F);
% 系统动力学方程
DynamicEqs = M*d2x+K*x==0;
% 第二问
% 求特征频率和特征向量
% 求特征值有两种方法，先对特征值做无量纲化，或者一直带量纲计算
% 如果不代入量纲一的参数，如下代码所示，有两个问题
% % 特征矩阵
% A = K-w^2*M;
% % 特征方程
% eqA = det(A);
% % 求解频率
% w_sol = solve(eqA,w)
% 一是解出来的频率是有可能带负号的，要区分一个符号表达式是否带负号还挺难的
% 使用量纲为一的参数，那么解得的频率的表达式都为正
% 二是解出来的频率可能为 RootOf 类型，需要使用 double 求得精确解
% 但是 double 不能对符号表达式使用，因此还需要先使用 subs 将频率的符号表达式中的 E I L m 代入一些数值
% 但是这样就失去了 E I L m，不够优雅
% 使用量纲为一的参数，那么解得的频率的表达式就不带 E I L m 了，就可以直接使用 double 了
% 其实 RootOf 类型的符号表达式也是可以被解决的
% 这样的话，理论上就可以用很少的代码求解出来，代码如下
% % 特征矩阵
% A = K-w^2*M;
% % 特征方程
% eqA = det(A);
% w_sol = solve(eqA,w,'MaxDegree', 3)
% % 特征向量矩阵初始化
% fai = zeros(3,size(w_sol,1));
% % 求特征向量
% for i = 1:1:size(w_sol,1)
% % 求解有量纲的特征向量_
% fai(:,i) = null(K-(w_sol(i))^2*M);
% end
% 但是这样有一个问题就是，对一元三次方程求出来的解，是带虚数的
% 这样求出来的特征值是没有意义的，但是我不知道怎么处理
% 所以也不能过分追求精确解
% 对代入了特征值的特征矩阵使用 null 函数，求得 empty，说明这个矩阵是满秩的
% 但是由于 |A-λE|=0，这个矩阵一定不是满秩的，说明前面的特征值求错了
% 立即可以想到 double 时一个近似的过程，使得类型转换后的特征值 λ1 与原特征值 λ 有偏差，使 |A-λ1E|=0 不成立了
% 想要在 double 之前判断 K-(w_sol(1))^2*M 的秩，会报错：
% 测试代码
% w_sol = (dPara*lambda_sol).^(1/2);
% class(K-(w_sol(1))^2*M)
% K-(w_sol(1))^2*M
% rank(K-(w_sol(1))^2*M)
% 错误信息
% 错误使用 symengine
% First argument must be a polynomial or polynomial expression.
% 出错 sym/privUnaryOp (第 1040 行)
% Csym = mupadmex(op,args{1}.s,varargin{:});
%
% 出错 sym/rank (第 13 行)
% r = double(privUnaryOp(A, 'symobj::rank'));
% 也就是说我想验证我解出来的原特征值 λ 是否能使 |A-λE|=0 成立都没办法
% 但是后面我灵机一动，突然想到，如果把 null 的参数转化成 double 类型的，再来试一次，会不会求出一个近似解
% 然后就可以求出来了，就很神奇
% 这说明，对于 sym 类型的矩阵，null 必须要矩阵不满秩才能求出来
% 但是如果不是 sym 类型的矩阵，null 的要求就会放松，允许求出一个近似解
% 综上，含无量纲化频率步骤的频率求解代码为
% 特征矩阵
A = K-w^2*M;
% 特征方程
eqA = det(A);
% 代入量纲一的参数 lambda = (m/k)*w^2
eqA = subs(eqA,w,(dPara*lambda)^(1/2));
% 求解无量纲的频率
% 解得的频率是从小到大排序的
lambda_sol = solve(eqA,lambda);
% 如果求得一个 RootOf 类型的符号式，说明需要使用 double 求精确解
try
% double 强制类型转换
lambda_sol = double(lambda_sol);
catch ME
% null
end
% 频率量纲化
w_sol = (dPara*lambda_sol).^(1/2);
% 特征向量矩阵初始化
fai = zeros(3,size(w_sol,1));
% 特征向量误差向量初始化
fai_error_norm = zeros(size(w_sol,1),1);
% 求特征向量
% 如果代入特征值的特征矩阵缺秩，可以直接使用 null
try
for i = 1:1:size(w_sol,1)
% 求解有量纲的特征向量
fai(:,i) = null(K-(w_sol(i))^2*M);
end
% 如果代入特征值的特征矩阵满秩，不能直接使用 null
catch ME
for i = 1:1:size(w_sol,1)
% 特征矩阵无量纲化，方便 null 求解
tmp = double((K-(w_sol(i))^2*M)/dPara/m);
% 求解特征向量
fai(:,i) = null(tmp);
% 计算特征向量的近似误差
fai_error_norm(i) = norm(tmp*fai(:,i));
end
end
% 第三问
% 近似频率初始化
w_prox_1 = 0;
% 邓克利法
for i = 1:1:size(F,1)
w_prox_1 = w_prox_1 + F(i,i)*M(i,i);
end
w_prox_1 = 1./sqrt(w_prox_1);
% 第四问
% 假设模态
fai_asum_2 = [[1 1.5 1]',[1 -1 -1]'];
% 近似频率初始化
w_prox_2 = sym(zeros(1,size(fai_asum_2,2)));
% 瑞利法
for i = 1:1:size(fai_asum_2,2)
w_prox_2(i) = sqrt(fai_asum_2(:,i)'*K*fai_asum_2(:,i)/(fai_asum_2(:,i)'*M*fai_asum_2(:,i)));
end
% 第五问
% 设一个频率，方便里兹法
syms w_3 lambda_3 real;
% 假设模态
fai_asum_3 = zeros(3,2,2);
fai_asum_3(:,:,1) = [1 1.5 1;-1.5 -2 -1.5]';
fai_asum_3(:,:,2) = [1 1.5 1;2 1 -1]';
% 近似频率初始化
w_prox_3 = sym(zeros(2,1,size(fai_asum_3,3)));
% 瑞利法
for i = 1:1:size(fai_asum_3,3)
% 缩并后的刚度矩阵
K_3 = fai_asum_3(:,:,i)'*K*fai_asum_3(:,:,i);
% 缩并后的质量矩阵
M_3 = fai_asum_3(:,:,i)'*M*fai_asum_3(:,:,i);
% 特征矩阵
A_3 = K_3-w_3^2*M_3;
% 特征方程
eqA_3 = det(A_3);
% 代入量纲一的参数 lambda = (m/k)*w^2
eqA_3 = subs(eqA_3,w_3,(dPara*lambda_3)^(1/2));
% 求解无量纲的频率
% 解得的频率是从小到大排序的
lambda_3_sol = solve(eqA_3,lambda_3);
% 如果求得一个 RootOf 类型的符号式，说明需要使用 double 求精确解
try
% double 强制类型转换
lambda_3_sol = double(lambda_3_sol);
catch ME
% null
end
% 频率量纲化
w_prox_3(:,1,i) = (dPara*lambda_3_sol).^(1/2);
end
% 第六问
% 设一个频率，方便里兹法
syms w_4 lambda_4 real;
% 假设模态
fai_asum_4 = zeros(3,2,2);
fai_asum_4(:,:,1) = [1 1.5 1;-1.5 -2 -1.5]';
fai_asum_4(:,:,2) = [1 1.5 1;2 1 -1]';
% 近似频率初始化
w_prox_4 = sym(zeros(2,1,size(fai_asum_4,3)));
% 瑞利法
for i = 1:1:size(fai_asum_4,3)
% 缩并后的动力矩阵
H_4 = fai_asum_4(:,:,i)'*M*F*M*fai_asum_4(:,:,i);
% 缩并后的质量矩阵
M_4 = fai_asum_4(:,:,i)'*M*fai_asum_4(:,:,i);
% 特征矩阵
A_4= M_4-w_4^2*H_4;
% 特征方程
eqA_4= det(A_4);
% 代入量纲一的参数 lambda = (m/k)*w^2
eqA_4= subs(eqA_4,w_4,(dPara*lambda_4)^(1/2));
% 求解无量纲的频率
% 解得的频率是从小到大排序的
lambda_4_sol = solve(eqA_4,lambda_4);
% 如果求得一个 RootOf 类型的符号式，说明需要使用 double 求精确解
try
% double 强制类型转换
lambda_4_sol = double(lambda_4_sol);
catch ME
% null
end
% 频率量纲化
w_prox_4(:,1,i) = (dPara*lambda_4_sol).^(1/2);
end
% 第七问
% 动力矩阵
D = F*M;
% 假设模态
fai_asum_5 = zeros(3,2,2);
fai_asum_5(:,:,1) = [1 1.5 1;-1.5 -2 -1.5]';
fai_asum_5(:,:,2) = [1 1.5 1;2 1 -1]';
% 矩阵迭代法
for i = 1:1:size(fai_asum_5,3)
% 一阶频率的 fai 初值
tmpfai1 = fai_asum_5(:,1,i);
% 初始化
A_old = zeros(size(tmpfai1,1),1);
A_new = zeros(size(tmpfai1,1),1);
res = inf;
% 迭代，直到误差小于给定精度
while res > 1e-7
% 如果没有初值，则赋初值为给定的 fai
if A_old == zeros(size(A_new,1),1)
% 先赋 A0
A_old = tmpfai1;
% 迭代前的 A 归一化
A_old = A_old/A_old(size(A_new,1));
end
% 迭代
A_new = D*A_old;
% 无量纲化
A_new = A_new*dPara;
% 迭代后的 A 归一化
A_new = A_new/A_new(size(A_new,1));
% 旧值更新
A_old = A_new;
% 计算误差
res = norm(A_new-A_old);
end
% 取 A_old 为第一阶振型
A_new = D*A_old;
% 使用 A_old 和 D*A_old 的最后一个元素计算基频
w_prox_5(1,1,i) = (A_old(size(A_old,1))/(A_new(size(A_old,1))))^(1/2);
% 求高阶振型
Mp1 = A_old'*M*A_old;
D1 = D-1/w_prox_5(1,1,i)^2/Mp1*(A_old*A_old')*M;
% 二阶频率的 fai 初值
tmpfai2 = fai_asum_5(:,2,i);
% 初始化
A_old = zeros(size(tmpfai2,1),1);
A_new = zeros(size(tmpfai2,1),1);
res = inf;
% 迭代，直到误差小于给定精度
while res > 1e-7
% 如果没有初值，则赋初值为给定的 fai
if A_old == zeros(size(A_new,1),1)
% 先赋 A0
A_old = tmpfai2;
% 迭代前的 A 归一化
A_old = A_old/A_old(size(A_new,1));
end
% 迭代
A_new = D1*A_old;
% 无量纲化
A_new = A_new*dPara;
% 迭代后的 A 归一化
A_new = A_new/A_new(size(A_new,1));
% 旧值更新
A_old = A_new;
% 计算误差
res = norm(A_new-A_old);
end
% 取 A_old 为第一阶振型
A_new = D1*A_old;
% 使用 A_old 和 D*A_old 的最后一个元素计算基频
% D1 不一定每项都是正数，这可能导致 A_old 或 A_new 的最后一项带负号
% 进而导致 w 的表达式中，出现根号下负号的情形 这使得 double 无法进行
% 因此需要先平方再开方消去负号
w_prox_5(2,1,i) = (A_old(size(A_old,1))/(A_new(size(A_old,1))))^2;
w_prox_5(2,1,i) = w_prox_5(2,1,i)^(1/4);
end
% 输出结果 第一问
% 柔度矩阵
disp('柔度矩阵');
disp(F);
% 刚度矩阵
disp('刚度矩阵');
disp(K);
% 质量矩阵
disp('质量矩阵');
disp(M);
% 系统动力学方程
disp('系统动力学方程');
disp(DynamicEqs);
% 输出结果 第二问
% 频率
disp('频率');
disp(w_sol);
% 简正模态
disp('简正模态');
disp(fai);
% 简正模态求解误差
disp('简正模态求解误差');
disp(fai_error_norm);
% 各阶简正模态形状
disp('各阶简正模态形状');
for i = 1:1:size(fai,1)
subplot(size(fai,1),1,i);
tmpx = 1:1:size(fai,1);
plot(tmpx,fai(:,i));
tmpstr = ['振型',num2str(i)];
title(tmpstr);
end
% 输出结果 第三问
% 邓克利法求系统基频
disp('邓克利法求系统基频为');
disp(w_prox_1);
% 输出结果 第四问
% 已知振型越接近真实振型，算出来的固有频率越准确
% 要衡量振型之间的接近程度，需要提供一套方法
% 一开始我以为是使用二范数，如下代码所示
% 已知 [1 1.8 2] 比 [1 1.4 1.8] 更接近真实振型，但是与真实振型的差的二范数反而更大
% fai_2 = [[0.4626 0.8608 1]',[-2.9339 -0.7458 1]',[1.4728 -3.1152 1]'];
% fai_asum = [1 1.4 1.8]';
%
% for i = 1:1:size(fai_2,2)
% fai_res = fai_asum - fai_2(:,i);
% res(1,i) = norm(fai_res);
% end
% double(res)
%
% fai_asum = [1 1.8 2]';
% for i = 1:1:size(fai_2,2)
% fai_res = fai_asum - fai_2(:,i);
% res(1,i) = norm(fai_res);
% end
% double(res)
%
% ans =
%
% 1.1043 4.5519 4.6098
%
%
% ans =
%
% 1.4734 4.7913 5.0381
% 于是使用第二种方法，计算向量之间的余弦值，如下代码所示
% fai_2 = [[0.4626 0.8608 1]',[-2.9339 -0.7458 1]',[1.4728 -3.1152 1]'];
% cosSim = @(a,b) sum(a.*b)/sqrt(sum(a.^2)*sum(b.^2));
%
% fai_asum = [1 1.4 1.8]';
% for i = 1:1:size(fai_2,2)
% cosRes(i) = cosSim(fai_asum, fai_2(:,i));
% end
% double(cosRes)
%
% fai_asum = [1 1.8 2]';
% for i = 1:1:size(fai_2,2)
% cosRes(i) = cosSim(fai_asum, fai_2(:,i));
% end
% double(cosRes)
%
% ans =
%
% 0.9960 -0.2744 -0.1218
%
%
% ans =
%
% 0.9996 -0.2487 -0.2073
% 这样算出来的接近度的可信度就很高了
% 瑞利法假设模态和对应的基频
for i = 1:1:size(w_prox_2,2)
disp('假设模态为');
disp(fai_asum_2(:,i));
disp('刚度矩阵表示的瑞利法求得对应的基频为');
disp(w_prox_2(i));
end
% 开始对比频率
disp('为了方便对比，下面将各频率无量纲化');
% 真实基频
w_1 = double(simplify(w_sol(1)/dPara^(1/2)));
disp('真实基频为');
disp(w_1)
% 瑞利法求得的无量纲基频初始化
w_2 = zeros(size(w_prox_2));
% 无量纲化
for i = 1:1:size(w_prox_2,2)
disp('刚度矩阵表示的瑞利法求得对应的基频为')
w_2(i) = double(simplify(w_prox_2(i)/dPara^(1/2)));
disp(w_2(i));
end
% 振型近似度
for i = 1:1:size(fai,2)
for j = 1:1:size(fai_asum_2,2)
cosRes(i,j) = cosSim(fai_asum_2(:,j), fai(:,i));
end
end
for j = 1:1:size(fai_asum_2,2)
str = ['第',num2str(j),'个假设模态相对真实模态的余弦近似度为'];
disp(str)
disp(double(cosRes(:,j)))
end
% 输出结果 第五问
% 里兹法假设模态和对应的频率
for i = 1:1:size(w_prox_3,3)
disp('假设模态为');
disp(fai_asum_3(:,:,i));
disp('刚度矩阵表示的里兹法求得对应的频率为');
disp(w_prox_3(:,i));
end
% 开始对比频率
disp('为了方便对比，下面将各频率无量纲化');
% 真实固有频率
for i = 1:1:2
tmpw = double(simplify(w_sol(i)/dPara^(1/2)));
tmpstr = ['第',num2str(i),'阶真实固有频率为'];
disp(tmpstr);
disp(tmpw)
end
% 里兹法求得的无量纲固有频率初始化
w_3 = zeros(size(w_prox_3,1),1,size(w_prox_3,3));
% 无量纲化
for i = 1:1:size(w_prox_3,3)
disp('刚度矩阵表示的里兹法求得对应的固有频率为');
w_3(:,1,i) = double(simplify(w_prox_3(:,1,i)/dPara^(1/2)));
disp(w_3(:,1,i));
end
% 振型近似度
for i = 1:1:size(fai_asum_3,3)
fai_tmp = fai_asum_3(:,:,i);
for j = 1:1:size(fai,2)
for j = 1:1:size(fai_asum_2,2)
cosRes(i,j) = cosSim(fai_tmp(:,j), fai(:,i));
end
end
for j = 1:1:size(fai_tmp,2)
str = ['第',num2str(j),'个假设模态相对真实模态的余弦近似度为'];
disp(str)
disp(double(cosRes(:,j)))
end
end
% 输出结果 第六问
% 里兹法假设模态和对应的频率
for i = 1:1:size(w_prox_4,3)
disp('假设模态为');
disp(fai_asum_4(:,:,i));
disp('柔度矩阵表示的里兹法求得对应的频率为');
disp(w_prox_4(:,i));
end
% 开始对比频率
disp('为了方便对比，下面将各频率无量纲化');
% 真实固有频率
for i = 1:1:2
tmpw = double(simplify(w_sol(i)/dPara^(1/2)));
tmpstr = ['第',num2str(i),'阶真实固有频率为'];
disp(tmpstr);
disp(tmpw);
end
% 里兹法求得的无量纲固有频率初始化
w_4 = zeros(size(w_prox_4,1),1,size(w_prox_4,3));
% 无量纲化
for i = 1:1:size(w_prox_4,3)
disp('柔度矩阵表示的里兹法求得对应的固有频率为');
w_4(:,1,i) = double(simplify(w_prox_4(:,1,i)/dPara^(1/2)));
disp(w_4(:,1,i));
end
% 输出结果 第七问
% 矩阵迭代法假设模态和对应的频率
for i = 1:1:size(w_prox_5,3)
disp('假设模态为');
disp(fai_asum_5(:,:,i));
disp('矩阵迭代法求得对应的频率为');
disp(w_prox_5(:,i));
end
% 开始对比频率
disp('为了方便对比，下面将各频率无量纲化');
% 真实固有频率
for i = 1:1:2
tmpw = double(simplify(w_sol(i)/dPara^(1/2)));
tmpstr = ['第',num2str(i),'阶真实固有频率为'];
disp(tmpstr);
disp(tmpw);
end
% 矩阵迭代法求得的无量纲固有频率初始化
w_5 = zeros(size(w_prox_5,1),1,size(w_prox_5,3));
% 无量纲化
for i = 1:1:size(w_prox_5,3)
disp('矩阵迭代法求得对应的固有频率为');
w_5(:,1,i) = double(simplify(w_prox_5(:,1,i)/dPara^(1/2)));
disp(w_5(:,1,i));
end
标签：disp,里兹法,end,asum,矩阵,fai,Matlab,迭代法,size 来源： https://blog.csdn.net/PriceCheap/article/details/122140338