我是靠谱客的博主 过时钢笔，最近开发中收集的这篇文章主要介绍[Matlab] 使用 solve 和 null 代替 eig 求解特征方程可能得到不一样的结果，觉得挺不错的，现在分享给大家，希望可以做个参考。
概述
1. 我以前写的啥我都忘了
比如在振动力学中求解阵型需要解特征方程
syms m1 m2 w l real;
>> syms g;
>> A=[-w^2*(3/2*m1+m2) -w^2*1/2*m2*l;-w^2*1/2*m2*l m2*g-w^2*3/4*m2*l^2];
>> solve(det(A),w)
警告: Solutions are only valid under certain conditions. To include parameters and conditions in the solution, specify the 'ReturnConditions' value as 'true'.
> 位置：sym/solve>warnIfParams (第 478 行)
位置: sym/solve (第 357 行)
ans =
0
-(4*((g*(3*m1 + 2*m2)*(9*m1 + 4*m2))/4)^(1/2))/(9*l*m1 + 4*l*m2)
(4*((g*(3*m1 + 2*m2)*(9*m1 + 4*m2))/4)^(1/2))/(9*l*m1 + 4*l*m2)
>>[V,D]=eig(A)
V =
[(- 3*l^2*w^2 + 4*g)/(2*l*w^2) + (2*((16*g^2*m2^2 - 24*g*l^2*m2^2*w^2 + 48*g*m1*m2*w^2 + 32*g*m2^2*w^2 + 9*l^4*m2^2*w^4 - 36*l^2*m1*m2*w^4 - 8*l^2*m2^2*w^4 + 36*m1^2*w^4 + 48*m1*m2*w^4 + 16*m2^2*w^4)^(1/2)/8 - (g*m2)/2 + (3*m1*w^2)/4 + (m2*w^2)/2 + (3*l^2*m2*w^2)/8))/(l*m2*w^2), (- 3*l^2*w^2 + 4*g)/(2*l*w^2) + (2*((3*m1*w^2)/4 - (16*g^2*m2^2 - 24*g*l^2*m2^2*w^2 + 48*g*m1*m2*w^2 + 32*g*m2^2*w^2 + 9*l^4*m2^2*w^4 - 36*l^2*m1*m2*w^4 - 8*l^2*m2^2*w^4 + 36*m1^2*w^4 + 48*m1*m2*w^4 + 16*m2^2*w^4)^(1/2)/8 - (g*m2)/2 + (m2*w^2)/2 + (3*l^2*m2*w^2)/8))/(l*m2*w^2)]
[ 1, 1]
D =
[(g*m2)/2 - (16*g^2*m2^2 - 24*g*l^2*m2^2*w^2 + 48*g*m1*m2*w^2 + 32*g*m2^2*w^2 + 9*l^4*m2^2*w^4 - 36*l^2*m1*m2*w^4 - 8*l^2*m2^2*w^4 + 36*m1^2*w^4 + 48*m1*m2*w^4 + 16*m2^2*w^4)^(1/2)/8 - (3*m1*w^2)/4 - (m2*w^2)/2 - (3*l^2*m2*w^2)/8, 0]
[ 0, (g*m2)/2 + (16*g^2*m2^2 - 24*g*l^2*m2^2*w^2 + 48*g*m1*m2*w^2 + 32*g*m2^2*w^2 + 9*l^4*m2^2*w^4 - 36*l^2*m1*m2*w^4 - 8*l^2*m2^2*w^4 + 36*m1^2*w^4 + 48*m1*m2*w^4 + 16*m2^2*w^4)^(1/2)/8 - (3*m1*w^2)/4 - (m2*w^2)/2 - (3*l^2*m2*w^2)/8]
>> simplify(D)
ans =
[(g*m2)/2 - (16*g^2*m2^2 - 24*g*l^2*m2^2*w^2 + 48*g*m1*m2*w^2 + 32*g*m2^2*w^2 + 9*l^4*m2^2*w^4 - 36*l^2*m1*m2*w^4 - 8*l^2*m2^2*w^4 + 36*m1^2*w^4 + 48*m1*m2*w^4 + 16*m2^2*w^4)^(1/2)/8 - (3*m1*w^2)/4 - (m2*w^2)/2 - (3*l^2*m2*w^2)/8, 0]
[ 0, (g*m2)/2 + (16*g^2*m2^2 - 24*g*l^2*m2^2*w^2 + 48*g*m1*m2*w^2 + 32*g*m2^2*w^2 + 9*l^4*m2^2*w^4 - 36*l^2*m1*m2*w^4 - 8*l^2*m2^2*w^4 + 36*m1^2*w^4 + 48*m1*m2*w^4 + 16*m2^2*w^4)^(1/2)/8 - (3*m1*w^2)/4 - (m2*w^2)/2 - (3*l^2*m2*w^2)/8]
>> simplify(V)
ans =
[(4*g*m2 + (16*g^2*m2^2 - 24*g*l^2*m2^2*w^2 + 48*g*m1*m2*w^2 + 32*g*m2^2*w^2 + 9*l^4*m2^2*w^4 - 36*l^2*m1*m2*w^4 - 8*l^2*m2^2*w^4 + 36*m1^2*w^4 + 48*m1*m2*w^4 + 16*m2^2*w^4)^(1/2) + 6*m1*w^2 + 4*m2*w^2 - 3*l^2*m2*w^2)/(4*l*m2*w^2), (4*g*m2 - (16*g^2*m2^2 - 24*g*l^2*m2^2*w^2 + 48*g*m1*m2*w^2 + 32*g*m2^2*w^2 + 9*l^4*m2^2*w^4 - 36*l^2*m1*m2*w^4 - 8*l^2*m2^2*w^4 + 36*m1^2*w^4 + 48*m1*m2*w^4 + 16*m2^2*w^4)^(1/2) + 6*m1*w^2 + 4*m2*w^2 - 3*l^2*m2*w^2)/(4*l*m2*w^2)]
[ 1, 1]
>>w1=(4*((g*(3*m1 + 2*m2)*(9*m1 + 4*m2))/4)^(1/2))/(9*l*m1 + 4*l*m2);
A1=[-w1^2*(3/2*m1+m2) -w1^2*1/2*m2*l;-w1^2*1/2*m2*l m2*g-w1^2*3/4*m2*l^2]
A1 =
[-(4*g*((3*m1)/2 + m2)*(3*m1 + 2*m2)*(9*m1 + 4*m2))/(9*l*m1 + 4*l*m2)^2, -(2*g*l*m2*(3*m1 + 2*m2)*(9*m1 + 4*m2))/(9*l*m1 + 4*l*m2)^2]
[ -(2*g*l*m2*(3*m1 + 2*m2)*(9*m1 + 4*m2))/(9*l*m1 + 4*l*m2)^2, g*m2 - (3*g*l^2*m2*(3*m1 + 2*m2)*(9*m1 + 4*m2))/(9*l*m1 + 4*l*m2)^2]
>> Z=null(A1)
Z =
-(l*m2)/(3*m1 + 2*m2)
1
可见，使用 solve 和 null 代替 eig，得到了一个形式简单的特征向量
噢……之后看了一下，其实不是 eig(K-w^2*M) 这么用的
应该是这样：
>> [V,D]=eig(K,M)
错误使用 sym/eig
输入参数太多。
这样 matlab 解不出来……这就很奇怪了
这样的话，有一个问题就是，不太能知道特征根是不是重根
比如
K =
[5*k, k, 2*k]
[ k, 3*k, 0]
[2*k, 0, 2*k]
>> M
M =
[2*m, -m, -m]
[ -m, 2*m, -m]
[ -m, -m, 2*m]
>> A2=K-w^2*M
A2 =
[- 2*m*w^2 + 5*k, m*w^2 + k, m*w^2 + 2*k]
[ m*w^2 + k, - 2*m*w^2 + 3*k, m*w^2]
[ m*w^2 + 2*k, m*w^2, - 2*m*w^2 + 2*k]
>> solve(det(A2),w)
ans =
(k*m)^(1/2)/m
-(k*m)^(1/2)/m
-(3^(1/2)*(k*m)^(1/2))/(3*m)
(3^(1/2)*(k*m)^(1/2))/(3*m)
这里，K 和 M 都是 3*3 的，但是只有两个特征根，那就是有一个重根，解题需要用到重根，但是我不知道哪个是重根
2. 重新实验
syms m1 m2 m3 k1 k2 k3 w real;
M=[m1 0 0;0 m2 0;0 0 m3];
K=[k1+k2 -k2 0;-k2 k2+k3 -k3;0 -k3 k3];
A = K-w^2*M;
% 法1
[V,D] = eig(A);
disp(V);
disp(D);
% 法2
w_sol = solve(det(A),w);
w_sol_positive = w_sol(size(w_sol,1)/2+1:size(w_sol,1));
for i = 1:1:size(w_sol_positive,1)
w_tmp = w_sol_positive(i);
A_tmp = K-w_tmp^2*M;
fai = null(A_tmp);
disp(fai);
end
disp(w_sol_positive);
测试出来起码 eig 还能有结果，solve 和 null 出来不一定有结果……
为了不费脑筋，还是 eig 吧……虽然 eig 出来的结果也不好看就是了
最后
以上就是过时钢笔为你收集整理的[Matlab] 使用 solve 和 null 代替 eig 求解特征方程可能得到不一样的结果的全部内容，希望文章能够帮你解决[Matlab] 使用 solve 和 null 代替 eig 求解特征方程可能得到不一样的结果所遇到的程序开发问题。
如果觉得靠谱客网站的内容还不错，欢迎将靠谱客网站推荐给程序员好友。
本图文内容来源于网友提供，作为学习参考使用，或来自网络收集整理，版权属于原作者所有。
![[Matlab] 使用 solve 和 null 代替 eig 求解特征方程可能得到不一样的结果](/uploads/reation/bcimg3.png)
发表评论 取消回复