标签：disp 里兹法 end asum 矩阵 fai Matlab 迭代法 size
《振动理论》课程考核
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
本站声明： 1. iCode9 技术分享网（下文简称本站）提供的所有内容，仅供技术学习、探讨和分享； 2. 关于本站的所有留言、评论、转载及引用，纯属内容发起人的个人观点，与本站观点和立场无关； 3. 关于本站的所有言论和文字，纯属内容发起人的个人观点，与本站观点和立场无关； 4. 本站文章均是网友提供，不完全保证技术分享内容的完整性、准确性、时效性、风险性和版权归属；如您发现该文章侵犯了您的权益，可联系我们第一时间进行删除； 5. 本站为非盈利性的个人网站，所有内容不会用来进行牟利，也不会利用任何形式的广告来间接获益，纯粹是为了广大技术爱好者提供技术内容和技术思想的分享性交流网站。