研究生课程论文封面
课程名称：数值分析
论文题目：基于Matlab对多自由度振动系统的数值分析学生班级；机械工程机械电子系3班
学生姓名：陈大爷
任课教师：王师傅
学位类别：学位课（2学分，32学时）
评分标准及分值选题与参阅资料
（分值）
论文内容
（分值）
论文表述
（分值）
创新性
（分值）
评分
论文评语：
总评分评阅教师: 评阅时间年月日
注：此表为每个学生的论文封面，请任课教师填写分项分值
基于Matlab 对多自由度振动系统的数值分析
摘要：多自由度主要研究矩阵的迭代求解，我们在分析抽象的理论的同时根据MATLAB 编程实现数据的迭代最后可以得到所要的数据，使我们的计算更加简便。

关键词： 振动系统；多自由度 ；迭代；MATLAB
引言：在工程振动中，研究某系统振动时，首先要求出系统的固有频率。

对于多自由度振动系统，计算系统固有频率与主振型主要有2种方法：(1)利用特征矩阵方程式与特征方程式求解；(2)矩阵迭代法求解【1】。

MATLAB 作为一个以矩阵和数组为核心计算的软件，对矩阵迭代法中的矩阵迭代计算尤其适合【2】。

本文主要利用MATLAB 对多自由度系统振动矩阵迭代求解。

一.多自由度振动系统 1.多自由度振动系统的数学模型
多自由度振动系统的数学模型【1】：
[]{}[]{}[]{}{}M x C x K x f ++= （1-1）
其中[]M 、[]C 、[]K 、{}f 和{}x 分别为质量矩阵、阻尼矩阵、刚度矩阵、力向量和响应向量。

把这个时域矩阵方程变换到拉氏域（变数为p ），并假定初始位移和初始速度为零，则得：
[][][]{}{}2()()()p M p C K X p F p ++= （1-2）
或 []{}{}()()()Z p X p F p = （1-3）
式中 []()Z p ：动刚度矩阵。

对于N 自由度系统，此方程有2N 个复共轭对出现的特征根：
i i i
i
i i j j λσωλσω*
=-+⎧⎨=--⎩ 其中i σ阻尼因子；i ω为阻尼固有频率。

将i φφλλ==和i 带入公式中，得
i i i φλφ=A （1-4）
二．多自由度振动系统的数值分析
1.多自由度振动系统的迭代运算
记1X 为初始迭代列阵，由展开定理，1X 可以表示为
111φa X =n n a a φφ+++ 22 （1-5）
对上式左乘矩阵A ，由式（1-4）得知第一次迭代后所得的列阵为
n n n a a a AX X φλφλφλ+++== 22211112=⎪⎪⎭⎫ ⎝
⎛+++n n n a a a φλλφλλφλ12122111 （1-6） 经第二次迭代后，得
⎥⎥⎦
⎤⎢⎢⎣
⎡⎪⎪⎭⎫ ⎝⎛++⎪⎪⎭⎫ ⎝⎛+==n n n a a a AX X φλλφλλφλ2
122
122112
123 同理第（r-1）次迭代后的结果为
⎥⎥⎦
⎤⎢⎢⎣⎡⎪⎪⎭⎫ ⎝⎛++⎪
⎪⎭⎫ ⎝⎛+==----n r n n r r r r a a a AX X φλλφλλφλ1
121
12
211111
（1-7）
2.矩阵迭代法的实例分析
图1-1 分析图
用矩阵迭代法求解过程如下：
解: 用影响系数法求得系统的质量矩阵和刚度矩阵为
⎥⎥
⎥⎦
⎤
⎢⎢⎢⎣⎡=m m m M 2000000⎥⎥⎥⎦
⎤
⎢⎢⎢⎣⎡----=K K
K K K
K K
K 220
2302 算出K 的逆阵及系统的动力矩阵为
⎥⎥⎥⎦
⎤
⎢⎢
⎢⎣⎡=-5.22122111111
k K ，⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡==-5214212111k m M K A 若[]T
X 1111=，第一次迭代后得到
⎥⎥⎥⎦
⎤
⎢⎢⎢⎣⎡==87411k
m AX Y ， ()⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡==000000.1857000.0500000.011312Y Y X
k
k
2k
m
2m
m
重复上述步骤，各次的迭代结果列于表1-1。

由表可见，经过6次迭代后已有87X X =，
所以第一阶主振型及基频取为⎥⎥⎥⎦
⎤
⎢⎢⎢⎣⎡==000000.1860806
.0462598.071X φ， m
k k
m 373087
.0184210
.711
1
1=
=
λω 表1-1 第一阶主振型的迭代
r
1
2
3
4
5
6
7
r X
1 1 1 0.500000 0.875000 1.000000 0.465517 0.862069 1.000000 0.462830 0.860911 1.000000 0.462617 0.860814 1.000000 0.462598 0.860806 1.000000 0.462598 0.860806 1.000000 1λm k 8
7.250000 7.189655 7.184652 7.184245 7.184210
三．利用MATLAB 对多自由度系统振动矩阵迭代求解
1.MATLAB 计算程序
clc; clear;
%建立质量矩阵M ，刚度矩阵K\ syms k m; M=[m 0 0 0 m 0 0 0 2*m]; K=[2*k -k 0 -k 3*k -2*k 0 -2*k 2*k];
%*********迭代第n 阶主阵型******************** %n 为计数器 % n=1;
%计算系统动力矩阵A
if n==1
A(:,:,n)=K\M;
elseif n~=1
A(:,:,n)=A(:,:,n-1)-(MP(:,n-1)\t(:,n-1))*f(:,n-1)*f(:,n-1)'*M;
end
%定义初始迭代向量X(1)
if n==1
X(:,1,n)=[1 1 1]';
elseif n==2
X(:,1,n)=[1 1 -1]';
elseif n==3
X(:,1,n)=[1 -1 1]';
end
%迭代过程，Y为中间矩阵,i为迭代次数
i=1;
Y(:,1,n)=A(:,:,n)*X(:,1,n);
X(:,2,n)=Y(3,1,n).\Y(:,1,n);
while
abs(X(1,i,n)-X(1,i+1,n))>=0.000001&&abs(X(2,i,n)-X(2,i+1,n))>0.000001 Y(:,i+1,n)=A(:,:,n)*X(:,i+1,n);
X(:,i+2,n)=Y(3,i+1,n).\Y(:,i+1,n);
X(:,2);
i=i+1;
end
% X(:,:,n)
f(:,n)=X(:,i,n);
t(n)=Y(3,i,n);
MP(:,n)=f(:,n)'*M*f(:,n);
end
%输出数据过程
disp('第一阶主阵型的迭代结果');
X(:,:,1)
disp('第二阶主阵型的迭代结果');
X(:,:,2)
disp('第三阶主阵型的迭代结果');
X(:,:,3)
disp('φi的计算结果，矩阵的每列分别是1,2 ，3阶的');
f
disp('######******注意***********主阵型的迭代结果后面的0是系统的占位符号，不算计算结果');
2.MATLAB编程计算结果
四．结束语
本文采用MATLAB语言，不仅相对语句少，可读性强，只要输入必要的数据，就可以快速地获得振动系统的固有频率以及主振型【3】。

对设计人员计算复杂多自由度系统固有频率具有参考意义。

参考文献
[1] 许本文,焦群英.机械振动与模态分析基础[M].北京:机械工业出版社,1988.
[2] 方治华,贾宏玉,沈利.应用Matlab对有死区的振动系统进行仿真分析[J],包头钢铁学院学报,2004,12.
[3 张志涌等,精通MA TLAB[M],北京:北京航空航天大学出版社,2000.。