[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]>> MM =[2*m, -m, -m]
[ -m, 2*m, -m]
[ -m, -m, 2*m]>> A2=K-w^2*MA2 =[- 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 求解特征方程可能得到不一样的结果相关推荐
- null在matlab中什么意思,null表示什么
Null 翻译为"空的".在计算机中常表示无结果,空值.是一种表示这个指标并不指向任何的对象的引用,它的元素只有"零",Null 是在计算中具有保留的值,用于指 ...
- 龙格库塔法解微分方程组的matlab程序,MATLAB实例源码教程：龙格库塔法求解微分方程组源代码实例.doc...
MATLAB实例源码教程:龙格库塔法求解微分方程组源代码实例.doc MATLAB实例源码教程龙格库塔法求解微分方程组源代码实例题目用经典 Runge-Kutta方法求下列一阶微分方程组的近似解y1 ...
- matlab用solve解方程错误提示,MATLAB中使用solve解决方程组的问题
希望使用MATLAB的solve函数解出一个带有虚数的方程组,但是一直提示计算错误,要么就是算不出来结果,希望大佬们能帮帮忙 程序如下: syms a1 a2; a=[a1 a2]; C11=3.06 ...
- help efun matlab,Matlab优化工具箱在函数最值求解中的应用.pdf
Matlab优化工具箱在函数最值求解中的应用.pdf 系 统 解 决 方 案 Matlab优化工具箱在函数最值求解中的应用 彭东海 (中山职业技术学院数学教研室,广 东 中山 528404) 摘 要 ...
- matlab怎么求三次微分,matlab课设三阶微分方程多种方法求解.doc
matlab课设三阶微分方程多种方法求解 目录 一.课程设计题目及意义 -------- 1 页 二.课程设计任务及要求 --------2 页 三.课程设计详细过程及结果 --------3至10页 ...
- 一阶欧拉近似matlab,MATLAB改进欧拉法与四阶龙格-库塔求解一阶常微分方程.doc
MATLAB改进欧拉法与四阶龙格-库塔求解一阶常微分方程 姓名:樊元君 学号:2012200902 日期:2012.11.06 一.实验目的 掌握改进欧拉法与四阶龙格-库塔求解一阶常微分方程的初值问题 ...
- matlab解方程大于0的解,matlab用solve解方程组，解出来有很多组解，如何编程只保留全部未知数都大于0的一组解（该方...
matlab用solve解方程组,解出来有很多组解,如何编程只保留全部未知数都大于0的一组解(该方 mip版 关注:130 答案:4 悬赏:40 解决时间 2021-01-25 20:34 已解 ...
- 用共轭梯度法求极小值matlab,用MATLAB实现最速下降法_牛顿法和共轭梯度法求解实例——张小强.doc...
机电产品优化设计 课程设计报告 姓 名 :张 小 强 学 号 :201222080633 学 院 :机械电子工程学院 实验的题目和要求 一．课程名称:最优化设计方法 二．实验日期:2013年6月27日 ...
- 【Matlab优化覆盖】虚拟力算法求解无线网络传感覆盖优化问题【含源码 1187期】
一.代码运行视频(哔哩哔哩) [Matlab优化覆盖]虚拟力算法求解无线网络传感覆盖优化问题[含源码 1187期] 二.matlab版本及参考文献 1 matlab版本 2014a 2 参考文献 [1 ...
最新文章
- 零起点学算法10——求圆柱体的表面积
- arcgis for android 学习 - (4) 了解mapView的一些方法和事件
- Qt实现应用程序单实例运行--LocalServer方式
- c++ template(4)基本技巧
- linux源码包编译安装与rpm安装方法介绍
- 从新手机到老股票 闲鱼为何会沦为骗子与营销的新平台？
- PRML-系列二之2.2
- 95-10-170-启动-KafkaRequestHandlerPool
- 计算机应用数学自考,计算机应用数学-补充题16年自考复习资料
- html设置等宽字体效果
- 世界各国国家代码简称
- sql 2008 R2 备份和还原
- Squid运行控制脚本_wuli大世界_新浪博客
- key位置 win10生成的ssh_WIN 10生成SSH密钥教程
- Community Preserving Network Embedding 论文笔记
- 隐藏通信隧道技术（下）
- 模板皮肤AnotherEon001中日历css自定义修改
- android照片视频备份,Android 保存图片或视频到相册并刷新相册
- Nico的刷题日记（三）
- 中级微观经济学：Chap 12 不确定性