理论上已经证明， 只要函数 f(t,y(t)) , y(tn ) 。
适当光滑，可以确认式（2-7）初值问题的解存在并且唯一的。 在本题中将采用 MATLAB 的 ode45 （） 函数对已经解耦的方程进行数值积分求解。 ode45 （）函数使用的是显式 Runge-Kutta(4,5)法则与 Dormand prince 组合解法的中阶解法， 所以特别适用于仿真线性化程度高的系统。由于 ode45（）函数计算快，一般来说，在解 题的第一次仿真时，采用 ode45（）函数是最佳选择。
2
哈尔滨工业大学课程设计
一 引言
机械动力学是机械原理的主要组成部分。它研究机械在运转过程中的受力、机械中 各构件的质量与机械运动之间的相互关系，是现代机械设计的理论基础。研究机械运转 过程中能量的平衡和分配关系。主要研究的是：在已知外力作用下，求具有确定惯性参 量的机械系统的真实运动规律、分析机械运动过程中各构件之间的相互作用力、研究回 转构件和机构平衡的理论和方法、机械振动的分析、以及机构的分析和综合等等。 由于本次作业主要研究的是“二自由度振动有阻尼及强迫振动响应”问题，所以本 文将先对该问题的两种解法进行相关的介绍，这两种方法分别是：解耦法和数值法。然 后根据课上给出的具体实例，应用这两种方法给出最终的结果。通过比较这两种结果， 加深我们对“二自由度振动有阻尼及强迫振动响应”问题的理解，并且在解决问题的过 程中，提高我们运用 MATLAB 的能力。
三 实例
3.1 问题描述
6
哈尔滨工业大学课程设计
3.2 解耦法求解
系统的图示为图（2程的矩阵表示：
9 0 2.7 0.3 27 3 0 0 1 x 0.3 0.3 x 3 3 x 1 sin 3t
式（3-3）
由式（3-3）可以看到，通过解耦法我们将原本耦合的二自由度有阻尼受迫振动系统 分解成了两个互不耦合的单自由度系统。利用坐标变换求出模态坐标下的初始条件后利 用 Matlab 里的 dsolve（）函数对微分方程进行求解，并用 simplify（）函数对求出的 解进行化简，我们可以得到模态坐标下的振动响应： r(1)= 0.0086*cos(3.0*t) + 0.1*sin(3.0*t) - (2.13*cos(1.41*t))/exp(0.1*t) (0.364*sin(1.41*t))/exp(0.1*t) r(2)= - 5.38*10^(-4)*cos(1.99*t)*(62.3*cos(1.01*t) - 2.65*cos(4.99*t) + 315.0*sin(1.01*t) - 66.1*sin(4.99*t)) - (2.09*cos(1.99*t))/exp(0.2*t) (0.00839*sin(1.99*t))/exp(0.2*t) - (sin(1.99*t)*(0.169*exp(0.2*t)*cos(1.01*t) + 0.0355*exp(0.2*t)*cos(4.99*t) - 0.0335*exp(0.2*t)*sin(1.01*t) 0.00142*exp(0.2*t)*sin(4.99*t)))/exp(0.2*t) 在得到模态坐标下的振动响应后，通过坐标变换将模态坐标下的振动响应变换到物
4．
化为：
P M 2 CM 2 P M 2 KM 2 P M

1

1

1

1

1 2
1 0 F1 0 1 F 2
式（2-5）
在方程式（2-5）两端同时左乘 PT ，方程化为：
PT M 2 CM 2 P PT M 2 KM 2 P PT M
式（2-1）
式（2-2）
将式（2-1）和式（2-2）写成矩阵形式，可以得到其振动的微分方程，用式（2-2） 表示。
1 0 F1 Mx Cx Kx BF (t ) 0 1 F2
式（2-3）
式中 M 为质量矩阵，C 为阻尼矩阵，K 为弹簧刚度矩阵。 解耦法的实质就是将微分方程式（2-3）的左侧矩阵都化成对角阵，这样就将两个耦 合的微分方程式（2-1）和式（2-2）化成两个互相不耦合的微分方程，大大简化了问题 的求解难度。在得出解耦的解后，还要利用坐标变换将解耦的解转换到物理坐标系下。 具体求解可分为八个步骤：
二 解法介绍
2.1 解耦法
图 2-1 二自由度振动有阻尼及强迫振动响应系统 如图 2-1 所示是一个典型的二自由度振动有阻尼及强迫振动响应系统， 其中 k1 、k 2 是 弹簧刚度系数， c1 、 c2 是阻尼系数， F1 、 F2 是两个质量块所受的外力， x1 、 x2 是两个质 量块的位移。 m1 、 m2 两个质量块在外力作用下做受迫振动，分别对两个质量块进行受力 分析，如图 2-2 所示。
2 做变量替换，令 x M q ，并在方程两侧同左乘 M 2 将微分方程的第一项系

1
1
1．
数 M 消掉。此时方程变为：
q M CM q M KM q M

1 2

1 2

1 2

1 2

1 2
1 0 F1 0 1 F 2
式（2-4）
4
哈尔滨工业大学课程设计
式（3-1）
在 Matlab 中使用函数 chol（）对质量矩阵 M 进行乔里斯基分解，从而求出 M 用关系式 x M q 对 x 进行变量替换，并且在方程两端同乘以 M 为：
1 2
1 2

1 2
。利
，则式（3-1）可以化
0.3 0.1 3 1 0 q q q sin 3t 0.1 0.3 1 3 1
2.2 数值积分法
dy (t ) f (t , y (t )) a t b y dt y (a) y0
式（2-7）
所谓数值积分解法就是通过数值积分的算法逐个求出区间[a,b]内若干个离散点
a t0 t1 tn b 处的近似值 y(t1 ), y(t2 ),
式（3-2）
为了继续对角化系数矩阵，我们可以先求出 q 的系数矩阵的特征向量 P，令 q Pr 做
变量替换。 然后方程两端同时左乘该特征向量的转置， 并同时右乘该特征向量， 则式 （3-2）
可以化为：
0
0.2
0 2 0 0.7071 sin 3t 0.4 0 4 0.7071
5． 6．
计算 S M P 和 S 的逆矩阵。 计算模态坐标下的初始条件 r (0) S 1 X 0 ， r (0) S 1 X 0 。

1 2
7．
在模态坐标下求解解耦后的振动响应 r (t ) 。
5
哈尔滨工业大学课程设计
8．
求解物理坐标下的振动相应 x(t ) Sr (t ) 。
7
哈尔滨工业大学课程设计
理坐标系下： x(1)= 1.27*10^(-4)*cos(1.99*t)*(62.3*cos(1.01*t) - 2.65*cos(4.99*t) + 315.0*sin(1.01*t) - 66.1*sin(4.99*t)) - 0.0236*sin(3.0*t) - 0.00203*cos(3.0*t) + (0.502*cos(1.41*t))/exp(0.1*t) + (0.0859*sin(1.41*t))/exp(0.1*t) + (0.492*cos(1.99*t))/exp(0.2*t) + (0.00198*sin(1.99*t))/exp(0.2*t) + (0.236*sin(1.99*t)*(0.169*exp(0.2*t)*cos(1.01*t) + 0.0355*exp(0.2*t)*cos(4.99*t) - 0.0335*exp(0.2*t)*sin(1.01*t) 0.00142*exp(0.2*t)*sin(4.99*t)))/exp(0.2*t) x(2)= (1.51*cos(1.41*t))/exp(0.1*t) - 0.0709*sin(3.0*t) 3.8*10^(-4)*cos(1.99*t)*(62.3*cos(1.01*t) - 2.65*cos(4.99*t) + 315.0*sin(1.01*t) - 66.1*sin(4.99*t)) - 0.00608*cos(3.0*t) + (0.258*sin(1.41*t))/exp(0.1*t) - (1.48*cos(1.99*t))/exp(0.2*t) (0.00593*sin(1.99*t))/exp(0.2*t) (0.707*sin(1.99*t)*(0.169*exp(0.2*t)*cos(1.01*t) + 0.0355*exp(0.2*t)*cos(4.99*t) - 0.0335*exp(0.2*t)*sin(1.01*t) 0.00142*exp(0.2*t)*sin(4.99*t)))/exp(0.2*t) 利用 Matlab 作图可得图（3-2） ，图中实线为图（2-1）中的质量块 1 的运动，图中 虚线为图（2-1）中的质量块 2 的运动。
3
哈尔滨工业大学课程设计
图 2-2 两个质量块受力分析 通过对两个质量块的受力分析，我们可以分别写出两个质量块振动的微分方程：
m1 x1 c1 x1 c2 ( x2 x1 ) k1 x1 k2 ( x2 x1 ) m2 x2 c2 ( x2 x1 ) k2 ( x2 x1 )
从式（2-4）可以看出微分方程的第一项系数已经被消掉。
2．
令 Kw M 2 KM

1

1 2
、 Cw M 2 CM

1