《振动理论》课程考核

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