第08课-Matlab机械振动建模分析方法-课件PPT

1、机械振动(Mechanical Vibration)交通与车辆工程学院交通与车辆工程学院 刚宪约刚宪约2007年年9月月27日日第二课 导论：Matlab基础2021/8/261主要内容 Matlab基本操作 特征值问题 积分变换 线性常微分方程主要内容 Matlab基本操作 特征值问题 积分变换 线性常微分方程Matlab基本操作 数组与矩阵操作 求解线性方程组 多项式运算 绘图操作 编程初步MATLAB的概况 MATLAB （MatrixLaboratory）是由美国MathWorks 公司开发的； MATLAB具有数值计算、符号计算和图形可视化三大功能； MATLAB的基本数据单位是矩阵

2、； 具有丰富的工具箱 源程序的开放性 程序解释执行机械振动(Mechanical Vibration)交通与车辆工程学院交通与车辆工程学院 刚宪约刚宪约2007年年9月月27日日MATLAB是一种高级程序语言Matrix Laboratory只需输入几行简单的命令，而不必用编程语言（C+、VB等）去编写一大串程序，就可以快速地得到结果并解决问题2021/8/266应用MATLAB解线性代数 Ax=B A,B,x是矩阵。当维数都是1时，则成为一般的线性方程。 A=3,B=6 x=B/A=2 A=5,12,8;6,5,8;9,6,10 B=7;11;7 x=AB=-4.7073;-1.2439;5

3、.6829矩阵的输入：行内用空格或逗号 行间用分号A=5 12 8; 6 5 8 ; 9,6,10基本矩阵ones(3), ones(2,5)zeros(3), zeros(2,5)eye(3), eye(3,4), eye(4,3)矩阵大小r,c=size(A), size(A)r、c分别为A的行数和列数命令的编辑修改方向键及控制键可以编辑修改已输入的命令 r1=(1+sqt(5)/2 sqrt命令超过一行，用连接下一行 s=1-1/2+1/3-1/4+1/5-1/6+1/7 . -1/8+1/9-1/10+1/11-1/12; 各符号间的空格可有可无，加上空格增加可读性。-1/12后的分号

4、使结果不回显常量和变量 缺省变量名ans pi i,j 复数 (-1)(1/2) Nan Not-A-Number 非数 0/0 inf 无穷大 1/0复数及其运算 表示：c1=1-2i, c2=3+sin(.5)*j 运算：c=c22 极坐标：b1=abs(c1), b2=angle(c1) 直角坐标：b3=real(c1), b4=imag(c1)c2 = 3.0000 + 0.4794ic = 8.7702 + 2.8766ib1 = 2.2361b2 = -1.1071（弧度）b21=b2*180/pi (角度）b3=1b4= -2矩阵运算 + - *A=1 2 3; 4 5 6; 7

5、 8 9 B=1 2 3; 4 5 6; 7 8 9C=A+B D=A-Bhelp fix fix(3.6) help rand rand(3,3) rand(3,3)A=fix(15*rand(4,3) B=fix(15*rand(4,4)C=A*B D=B*A矩阵运算 invA=fix(10+17i) *rand(3,3) E=AA=5,12,8;6,5,8;9,6,10 B=7;11;7x=AB=-4.7073;-1.2439;5.6829x=inv(A)*B简单数组的访问形式 x=0 0.1*pi 0.2*pi 0.3*pi 0.4*pi 0.5*pix = 0 0.3142 0.62

6、83 0.9425 1.2566 1.5708 y=sin(x)y = 0 0.3090 0.5878 0.8090 0.9511 1.0000 x(3)ans = 0.6283x= 0 0.3142 0.6283 0.9425 1.2566 1.5708y= 0 0.3090 0.5878 0.8090 0.9511 1.0000 x(1:5)ans = 0 0.3142 0.6283 0.9425 1.2566 x(2:2:7)ans = 0.3142 0.9425 1.5708 y(5 2 6 1)ans = 0.9511 0.3090 1.0000 0如何创建数组 X=a b c d;

7、 由指定元素建立行向量 x=初值：终值 从初值到终值，增量为1 x=初值：增量：终值 x=linspace(初值,终值,n) 指定元素个数，线性等距分布 x=linspace(1,100,10) x=linspace(0,100,10) x=linspace(0,100,9) x=linspace(0,100,11) x=logspace(初值,终值,n) x=logspace(0,2,11)指定元素个数，对数等距分布数组的格式 a=1,2,3,4,5: 行向量 a=1 2 3 4 5: 行向量 b=1;2;3;4;5: 列向量 c=1 2 3;4 5 6;7 8 9: 矩阵 d=c 共轭转置

8、纯量与数组的运算 d=1 2 3;4 5 6;7 8 9; 2*d-1 ans = 1 3 5 7 9 11 13 15 17数组之间的运算 b=1 2 3 c=1;2;3b*c 与b.*c的区别b*c = 14c*b=1 2 3; 2 4 6; 3 6 9b.*c? Error using = .*Matrix dimensions must agree.c.*b? Error using = .*Matrix dimensions must agree.b.*c = 1 4 9c.*b = 1; 4; 9数组操作 A(r,c)：指定某一元素 A(r,:)：指定某一行元素 A(:,c)：指定

9、某一列元素A=1 2 3;4 5 6;7 8 9;A(1,2)A(2,:)A(:,3)数组规模 r,c=size(A)r、c分别为A的行数和列数A=1 2 3;4 5 6;7 8 9 r,c=size(A) A(2,:) n=length(b)n为向量b的长度 n=length(A(2,:) 多项式求根 多项式表示法 用行向量P表示 求已知多项式的根 r=roots(P) 求一组根所对应的多项式 P=ploy(r) P=1 -12 0 25 116;表示多项式 x4-12x3+25x+116P=1 -12 0 25 116;r=roots(P)r = 11.7473 2.7028 -1.225

10、1 + 1.4672i -1.2251 - 1.4672ip1=poly(r)p1 = 1.0000 -12.0000 0 25.0000 116.0000多项式的基本运算 多项式相乘c=conv(a,b) 如: a(x)=x3+2x2+3x+4 b(x)= x3+4x2+9x+16求c(x)=a(x)*b(x)= x6+6x5+20 x4+50 x3+75x2+84x+64 多项式相加、减，即矩阵相加，但规模需相同如a(x)=x3+2x2+x+1;b(x)= x+4求c(x)=a(x)+b(x) a=1 2 3 4; b=1 4 9 16; c=conv(a,b) a=1 2 1 1; b=

11、0 0 1 4; c=a+b;多项式的基本运算 多项式项除 q,r=deconv(a,b) 两多项式项除，q为商、r为余项 多项式微分 ployder(p)44*0 72*1 69*248*3 20*4 6*5 1*6 q,r=deconv(1 2 1,1 1) q = 1 1 r = 0 0 0 polyder(1 6 20 48 69 72 44) ans = 6 30 80 144 138 72基本二维画图命令Plot的用法 plot(x,y) plot(y) plot(x1,y1,x2,y2)x为横轴，y为纵轴x=linspace(0,2*pi,30);y=sin(x);z=cos(x

12、);plot(y);pause;plot(x,y);pause;plot(x,y,x,z)01234567-1-0.8-0.6-0.4-0.200.20.40.60.81图形标注 Title(text)：图形顶部加标题 xlabel (text)：横轴标注 ylabel (text)：纵轴标注 text(x,y,text)： (x,y)处加标注 gtext(text): 鼠标点击处加标注图形标注举例x=linspace(0,2*pi,30);y=sin(x);plot(x,y)title(y=sin(x)xlabel(x)ylabel(sin(x)text(pi/2,0.8, max)gtex

13、t(min)01234567-1-0.8-0.6-0.4-0.200.20.40.60.81y=sin(x)xsin(x)maxmin图例说明 Legend(s1,s2)建立一个图例说明盒在右上角 legend offplot(x,y,b,x,z,r)legend(sin(x),cos(x)pauselegend off01234567-1-0.8-0.6-0.4-0.200.20.40.60.81sin(x)cos(x)其它二维图形命令 bar条形图 例：x=-2.9:0.2:2.9;y=exp(-x.*x);bar(x,y)-3-2-1012300.10.20.30.40.50.60.70

14、.80.91其它二维图形命令 stairs 阶梯图 例：x=-2.9:0.2:2.9;y=exp(-x.*x);stairs(x,y)-3-2-1012300.10.20.30.40.50.60.70.80.91其它二维图形命令 stem火柴杆图 例：x=linspace(0,4*pi);y=exp(-0.3*x).*sin(x);stem(x,y);title(stem plot)02468101214-0.3-0.2-0.100.10.20.30.40.50.60.7stem plot基本三维曲线命令 plot3(x,y,z)x,y,z为同规模的向量或矩阵 plot3(x,y,x,s)字符

15、串s指定线型、颜色、点标记 plot3(x1,y1,z1,s1,x2,y2,z2,s2,)基本三维曲线命令举例 绘制三维圆柱螺线t=linspace(0,10*pi);x=sin(t)y=cos(t)z=tplot3(sin(t),cos(t),t)-1-0.500.51-1-0.500.51010203040曲面彩色网线图 mesh(x,y,z)x(n维)向量,y(m维)向量,z(mXn维)矩阵绘制网线图 例：x=-7.5:0.5:7.5; y=x;x,y=meshgrid(x,y);r=sqrt(x.2+y.2)+eps;z=sin(r)./r;mesh(x,y,z)2222/ )sin(

16、yxyxz曲面彩色曲面图x=-7.5:0.5:7.5; y=x;x,y=meshgrid(x,y);r=sqrt(x.2+y.2)+eps;z=sin(r)./r;mesh(x,y,z)figuresurf(x,y,z)2.3 MATLAB 语言流程控制 循环结构 for 结构 while 结构【例2-12】用循环求解【例2-13】用循环求解 求最小的 m【例2-14】求转移结构【例2-15】用循环求解 求最大的 m2.3.3 开关结构2.4.1 MATLAB 语言的函数的基本结构 nargin， nargout, varargin, varargout一个有趣的问题 房间里面有序号分别为1，

17、2，100的电灯，初始时都关闭； 又有100个序号分别为1，2，100的学生，当他们依次从这100盏灯前走过时分别将序号是本人序号的整数倍的所有灯开关切换一次（原来为开的变为关，原来为关的变为开）。 请问当所有的人都走过之后，哪些灯是亮的？主要内容 Matlab基本操作 特征值问题 积分变换 线性常微分方程特征值问题(1) 线性代数上所学的特征值问题 广义特征值问题xAxBxAx d = eig(A) returns a vector of the eigenvalues of matrix A. d = eig(A,B) returns a vector containing the gen

18、eralized eigenvalues, if A and B are square matrices. V,D = eig(A) produces matrices of eigenvalues (D) and eigenvectors (V) of matrix A, so that A*V = V*D. V,D = eig(A,B) produces a diagonal matrix D of eneralized eigenvalues and a full matrix V whosecolumns are the corresponding eigenvectors so th

19、at A*V = B*V*D . A = 1.6294 0.6324 0.9575 0.9572 0.6324 0.1951 0.9649 0.4854 0.9575 0.9649 0.3152 0.8003 0.9572 0.4854 0.8003 0.2838 V D=eig(A) V = -0.0399 -0.2694 0.6711 0.6896 -0.5964 -0.4213 -0.5784 0.3637 0.7572 -0.0454 -0.4537 0.4676 -0.2632 0.8648 -0.0966 0.4166 D = -0.7734 0 0 0 0 -0.2929 0 0

20、 0 0 0.2993 0 0 0 0 3.1905主要内容 Matlab基本操作 特征值问题 积分变换 线性常微分方程积分变换 Fourier变换 离散Fourier变换Fourier生平生平 1768年生于法国年生于法国 1807年提出年提出“任何周任何周期信号都可用正弦函期信号都可用正弦函数的级数表示数的级数表示” 1822年发表年发表“热的分热的分析理论析理论”，首次提出，首次提出“任何非周期信号都任何非周期信号都可用正弦函数的积分可用正弦函数的积分表示表示”Fourier变换(1) Fourier series: Any periodic motion can be represen

21、ted by a series of sines and cosines that are harmonically related.Fourier变换(2)10sincos2)(nnnnntbtaatx2/2/cos)(2TTnntdttxTa2/2/sin)(2TTnntdttxTb2/2/0)(2TTdttxTa正确吗？正确吗？Fourier变换(3) Exponential form of Fourier series:ntjnnntjnntjnnntjnnntjnnnecececaejbaejbaatx1*0102 21212)(2/2/)(1TTtjnndtetxTcFourier

22、变换(4) 如果函数的周期趋向于无穷大，对于Fourier级数应该做何变化？ 可以令TcXnnntjnnectx)(2/2/)(1TTtjnndtetxTcdeXtxeXeXTtxtjntjnnntjnn )(21)(21)(05.2.2 Fourier 变换的计算机求解傅立叶变换 傅立叶变换的意义 数学意义 从一个函数空间(集合)到另一个函数空间(集合)的映射； f(x)称为变换的原函数(相当于自变量)，F()称为象函数。 应用意义 把任意函数分解为简单周期函数之和，F()的自变量为频率，函数值为对应的振幅。 物理意义 把一般运动分解为简谐运动的叠加； 把一般电磁波(光)分解为单色电磁波(光

23、)的叠加。离散傅里叶变换 离散傅里叶变换是信号分析中的一种重要工具，它将时域内的问题转化为频域内的问题，在很多情况下大大地简化了问题的求解过程；另外计算时域信号的频谱也主要依靠离散傅里叶变换来完成。离散序列x（k）的离散傅里叶变换（DFT）定义为：1, 2 , 1 , 0)()(10)/2(NnekxnxNknkNj离散序列x（n）的离散傅里叶逆变换定义为：1, 2 , 1 , 0)(1)(10)/2(NkenxNkxNnnkNj快速傅立叶变换 快速傅立叶变换（FFT）是离散傅里叶变换的一种快速算法。（1）fft函数：一维快速傅立叶变换 调用形式：Y=fft（X，N） 返回N点的离散傅里叶变换

24、（2）ifft函数：一维快速傅立叶逆变换 调用形式：Y=ifft（X，N） 注：（1）当N取2的整数幂时，傅立叶变换的计算速度最快； 通常取大于且最靠近X实际长度的幂次。 N =2p p=nextpow2（X的实际长度）（2）一般情况下，fft的结果为复数，可用abs及angle 分别求其幅度、相位。例1、已知模拟信号)8cos(7)6sin(2)(tttx求N=128点DFT的幅度谱和相位谱。MATLAB命令窗口输入如下：N=128; n=0:127;t=0:0.01:127; %时域数据点数q=n*2*pi/N; %频域数据点数X=2*sin(6*pi*t)+7*cos(8*pi*t);Y

25、=fft(X,N);subplot(2,1,1);plot(q,abs(Y)subplot(2,1,2);plot(q,angle(Y)主要内容 Matlab基本操作 特征值问题 积分变换 线性常微分方程线性常微分方程(1)微分方程的解析解方法微分方程的解析解方法线性微分方程的数值解Matlab函数 T,Y = solver(odefun,tspan,y0) T,Y =solver(odefun,tspan,y0,options) where solver is one of ode45, ode23, ode113, ode15s, ode23s, ode23t,or ode23tb TCo

26、lumn vector of time points YSolution array. Each row in y corresponds to the solution at a time returned in the corresponding row of t.算例txxx5 . 0sin10010010 function xdot=testf(t,x)xdot(1,1)=x(2);xdot(2,1)=-x(2)/10-10*x(1)+sin(pi*t/2);TOUT,YOUT = ode45(testf,0 10,0 0);plot(TOUT,YOUT(:,1)MATLAB常用的基本数学函数常用的基本数学函数abs(x)：纯量的绝对值或向量的长度 angle(z)：复 数z的相角(Phase angle) sqrt(x)：开平方 real(z)：复数z的实部 imag(z)：复数z的虚 部 conj(z)：复数z的共轭复数 round(x)：四舍五入至最近整数 fix(x)：无论正负，舍去小数至最近整数 floor(x)：地板函数，即舍去正小数至最近整数 ceil(x)：天花板函数，即加入正小数至最近整数 rat(x)：将实数x化为分数表示 rats(x)：将实数x化为