建立系统的动力学方程
建立系统的动力学方程的方法：1、牛顿力学：牛顿第二定律；2、分析力学：拉格朗日方程。 以x、θ为两个变量建立二自由度系统动力学方程；3、影响系数法-张量算法 1、 根据牛顿第二定律，
2、 拉格朗日方程，
d∂L idt∂q
-∂L∂qi
=Qi(i=1,2,.....n) -------------拉格朗日第二类方程
L=T-V，称为拉氏函数，泛函；
势能函数：V=V(q1,q2,.......qn)；
1,q 2,.......q n)； 动能函数：T=T(q
广义坐标qi(i=1,2,.....n)对应的非保守力：Qi(i=1,2,.....n)
d∂L idt∂q
-∂L∂qi
=0(i=1,2,.....n) -------------保守系统的拉氏方程
-------------------利用上诉拉氏方程求解-------------------------------------
T=T=V=
121212
TMq ；V=q(mx
2
2
+Jθ )
12
qKq
T
[k1(x-l1θ)+k2(x+l2θ)]
22
将L=T-V代入拉氏方程可解的：
可见：与牛顿第二定律求得的结果一致。
+Kq=Q;对保守系统Q=0； M q
n
∑(m
j=1
ij
j+kijqj)=Qi(i=1,2,.....n) q
M=(
m0
0J
)
K=(
k1+k2k2l2-k1l1
k2l2-k1l1k1l1+k2l2
2
2
)
q=() θ
j+kijqj)=Qi(i=1,2,.....n)有明确的物理意义：弹性恢复力 +Kq=Q 和∑(mijqM q
j=1n
x
与保守力Q平衡。如张量理论，mij、kij可认为是张量的坐标，kij表-Kq、惯性力-M q
示：使系统仅产生沿qj坐标的单位位移时，沿qi坐标必须施加的外力Qi，或者说Q的i分量在q的j分量上的影响量的投影。
3、影响系数法-张量算法
！！！！！！！下面由张量分量投影计算理论直接求解-----------！！！！！！！！ 设
m11 * + k11 k12 * x
m21 θ + k21 k22 θ
k11:仅当x动1单位时，x向的作用力为k1*1+k2*1=k1+k2;
k12:仅当θ动1单位时，x向的作用力为-k1*(l1*1)+k2*(l2*1)=-k1*l1+k2*l2;
k21:仅当x动1单位时，θ向的作用力为-(k1*1)*l1+(k2*1)*l2=-k1*l1+k2*l2;
k22:仅当θ动1单位时，θ向的作用力为 k1*(l1*1)*l1+k2*(l2*1)*l2=k1*l1^2+k2*l2^2;
m11:仅当x”动1单位时，x”向的作用力为m*1;
k12:仅当θ”动1单位时，x”向的作用力为0; %仅绕质心转动时不影响x”向惯性力 k21:仅当x”动1单位时，θ”向的作用力为0; %仅平动时不影响θ”向惯性力？？？？？？不过质心时该怎么计算？？若旋转中心偏离质心a,则变为
ma
此时，k12: 仅当θ”动1单位时，x”向的作用力为m*(a*1); k21:仅当x”动1单位时，θ”向的作用力为m*1*a;
k22:仅当θ”动1单位时，θ”向的作用力为 (J+m*a^2)*1;
k22:仅当θ”动1单位时，θ”向的作用力为 J*1;
可见与上述结果一致。
对建立的动力学方程更换坐标求偏频
对上述系统建立前后轮纵向位移x1、x2的动力学方程 :x1=x-l1*θ；x2=x+l2*θ。使用matlab的solve（‘x1=x-l1*θ；x2=x+l2*θ’，‘x1’, ‘x2’）可以直接解出: θ=(x1-x2)/(l1+l2);x=(x1*l2+l1*x2)/(l1+l2),代入前面创建的方程组：
消去x1和x2,可得到如下的方程：
ml2+ρl
2
22
2
1+mx
l1l2-ρl
2
2
2+k1x1=0x
m
l1+ρl
2
2
2+mx
l1l2-ρl
2
2
1+k2x2=0x
⎡l22+ρ2⎢m2
l
所以⎢2
⎢ml1l2-ρ
2⎢l⎣
2
2
l1l2-ρ⎤m⎥⎡ 2 1⎤⎡k1xl⨯⎥⎢⎥+⎢22
2⎦⎣0xl+ρ⎥⎣
m12
⎥l⎦
0⎤⎡x1⎤
⎥⨯⎢⎥=0 k2⎦⎣x2⎦
1+η1* 2+ω1*x1=0xx
2
2+η2* 1+ω2*x2=0xx
l1l2-ρ
2
式中，η1=l2+ρ2 联系系数，表示两坐标之间的联系
2η2=
l1l2-ρl1+ρ
2
22
2
2
ω1=
k1l
2
m(l2+ρ)
偏频，表示前后悬挂独立振动时的振动频率，即x1=0时的振动频
率是w2,x2=0时的振动频率是w1,不同于系统的固有频率（2自由度独立时才相等）。 ω2=
22
m(l1+ρ)
k2l
2
ρ=
J
m 汽车绕质心轴的回转半径
在汽车设计中，希望行车时一个悬挂的振动不传到另一个悬挂上，为此，应使车身质量分布和前后轮位置满足：
质量分配系数ε=ll=1，这时η1=η2=0
12ω1=
k1lml2k2lml1
2
ρ
2
=
m1 k2
k1
ω2=
=
m2
2
2
J=mρ
=m1l1+m2l2
ε≠1
对于一般质量分配系数的耦合情况，可以用模态分析法求固有频率及其通解：
x1=A11sin(p1t+ϕ1)+A12sin(p2t+ϕ2)x2=A21sin(p1t+ϕ1)+A22sin(p2t+ϕ2)=β1A11sin(p1t+ϕ1)+β2A12sin(p2t+ϕ2)
⎡x1⎤⎡1即⎢x⎥=⎢β⎣2⎦⎣11⎤⎡A11sin(p1t+ϕ1)⎤⎡1⎥⎢⎥=⎢β2⎦⎣A12sin(p2t+ϕ2)⎦⎣β11⎤⎡xp1⎤
⎥⎥⎢
β2⎦⎣xp2⎦
可见，特征向量阵(模态矩阵)ф组成坐标变换矩阵（由老基到新的主坐标基的坐标变换矩阵），xp=фx为主坐标，主振动的坐标，在该坐标系上，各自由度独立振动（解耦）。
⎡1⎤⎡1⎤
⎢⎥⎢⎥⎣β1⎦⎣β2⎦
T
分别是新的主坐标基的两个基矢量在老基下的投影坐标。
在新的主坐标基下，Mp=фMф为主质量阵（主质量组成的对角阵）；Kp=фKф为主刚度阵（主刚度组成的对角阵）。这种矩阵变换的本质是张量的坐标变换。
+Kpx=0 主坐标方程组为解耦方程组。 Mp x
T
T
利用特征值分解找到系统的主坐标基，通过坐标变换进行解耦、简化计算、然后再变换回去，这是坐标变换的意义所在。
用matlab特征值分解法求平等与转动主模态（振型）
%SH760小轿车空载主要参数 m=1340; a=1.54;
b=1.29;
Ic=2395; %绕质心的转动惯量 k1=40000; k2=44000; M=[m,0;0,Ic];
K=[k1+k2,-(k1*a-k2*b);-(k1*a-k2*b),k1*a^2+k2*b^2];
[eig_vec,eig_val] = eig(inv(M)*K);
[omeg,w_order] = sort(sqrt(diag(eig_val))); %频率 mode_vec = eig_vec(:,w_order); %振型 T=2.*pi./omeg; %周期
mode_vec(:,1)=mode_vec(:,1)./mode_vec(1,1); mode_vec(:,2)=mode_vec(:,2)./mode_vec(1,2); subplot(2,1,1)
plot([1;2],mode_vec(:,1))
title(strcat('w1=',num2str(omeg(1)))); subplot(2,1,2)
plot([1;2],mode_vec(:,2))
title(strcat('w2=',num2str(omeg(2))));
因为对特征值进行了排序，所以w1
求平动与平动主模态（振型）与解析法仿真计算
%SH760小轿车空载主要参数 clear; m=1340; a=1.54; b=1.29;
l=a+b;
Ic=2395; %绕质心的转动惯量 rou=sqrt(Ic/m); k1=40*1000; k2=44*1000;
M=[m*(b^2+rou^2)/l^2,m*(a*b-rou^2)/l^2;m*(a*b-rou^2)/l^2,m*(a^2+rou^2)/l^2]; K=[k1,0;0,k2];
%用matlab特征值分解法求主振型------------------------------------------------
[eig_vec,eig_val] = eig(inv(M)*K);
[omeg] = (sqrt(diag(eig_val))); %频率不用sort排序 mode_vec = eig_vec;%(:,w_order); %振型 T=2.*pi./omeg; %周期
mode_vec(:,1)=mode_vec(:,1)./mode_vec(1,1); mode_vec(:,2)=mode_vec(:,2)./mode_vec(1,2);
w1=sqrt((k1*l^2)/(m*(b^2+rou^2))); w2=sqrt((k2*l^2)/(m*(a^2+rou^2))); w1_pian=sqrt((k1*l)/(m*b)); w2_pian=sqrt((k2*l)/(m*a));
subplot(2,2,1)
plot([1;2],mode_vec(:,1))
title(strcat('w_1=',num2str(omeg(1)),';w_1pian=',num2str(w1_pian))); subplot(2,2,2)
plot([1;2],mode_vec(:,2))
title(strcat('w_2=',num2str(omeg(2)),';w_2pian=',num2str(w2_pian)));
%完全用matlab sym解析法运算求解------------------------------------------------
syms A11 A12 phi1 phi2 t
x1=A11*sin(omeg(1)*t+phi1)+A12*sin(omeg(2)*t+phi2);
x2=mode_vec(2,1)*A11*sin(omeg(1)*t+phi1)+mode_vec(2,2)*A12*sin(omeg(2)*t+phi2); dx1=diff(x1); dx2=diff(x2);
x1_0=subs(x1,'t',0); x2_0=subs(x2,'t',0); dx1_0=subs(dx1,'t',0);
dx2_0=subs(dx2,'t',0);
eq=[sym(strcat(char(x1_0),'=1'));sym(strcat(char(dx1_0),'=0')); sym(strcat(char(x2_0),'=0'));sym(strcat(char(dx2_0),'=0'))];
s=solve_sym(eq);
x1=s.A11(1)*sin(omeg(1)*t+s.phi1(1))+s.A12(1)*sin(omeg(2)*t+s.phi2(1));
x2=mode_vec(2,1)*s.A11(1)*sin(omeg(1)*t+s.phi1(1))+mode_vec(2,2)*s.A12(1)*sin(omeg(2)*t+s.phi2(1)); ti=0:0.02:10;
x1i=subs(x1,'t',ti); x2i=subs(x2,'t',ti);
subplot(2,2,3)
plot(ti',[x1i',x2i'])
%用matlab指数运算求解----------------------------------------------------------
x0=[1;0];xd0=[0;0]; %初始条件 tf=10;dt=0.02; %时间向量
A=[zeros(2,2),eye(2);-M\K,zeros(2,2)]; Y=[x1;x2;x1';x2']
%expm(A)的意义是将坐标先变换到主坐标系，对对角值进行exp运算后再变换到原坐标系，如同张量坐标变换help expm
y0=[x0;xd0]; %四元变量的初始条件
for i=1:round(tf/dt)+1 %设定计算点，作循环计算 tj(i)=dt*(i-1);
y(:,i)=expm(A*tj(i))*y0; %循环计算矩阵指数 end
subplot(2,2,4),plot(tj,[y(1,:)',y(2,:)']),grid
%四阶参数矩阵
Y'=AY-->Y=expm(A*t)*Y0
可见，坐标的选取对固有频率没有影响，但对振型有影响。W1_pianpin和w2_pianpin是前后的偏频（假设质量分配系数为1计算）。
用matlab ode45()直接进行仿真计算
%SH760小轿车空载主要参数 clear; m=1340; a=1.54; b=1.29; l=a+b;
Ic=2395; %绕质心的转动惯量 rou=sqrt(Ic/m); k1=40*1000; k2=44*1000;
M=[m*(b^2+rou^2)/l^2,m*(a*b-rou^2)/l^2;m*(a*b-rou^2)/l^2,m*(a^2+rou^2)/l^2]; K=[k1,0;0,k2];
%用matlab ode45数值解------------------------------------------------
A=[zeros(2,2),eye(2);-M\K,zeros(2,2)]; %四阶参数4X4矩阵X'=AX-->X=expm(A*t)*X0 X=[x1;x2;x1';x2']
syms x1 x2 dx1 dx2
df_sym=A*[x1;x2;dx1;dx2];
df_sym=subs(df_sym,{'x1','x2','dx1','dx2'},{'x(1)','x(2)','x(3)','x(4)'}); n=length(df_sym);
i=1;
ss='['; %先定义好很重要，否则再循环体中定义时，每一循环ss不累加。 while i
ss=strcat(ss,char(df_sym(i)),';'); i=i+1;
end
ss=strcat(ss,char(df_sym(i))); ss=strcat(ss,']'); f=inline(ss,'t','x');
[t,x]=ode45(f,[0 10],[1,0,0,0]);%初始y=0,y'=1 %subplot(2,2,1)
plot(t,[x(:,1),x(:,2)]) %时间状态系列
用s-function进行仿真计算
%sh760.m
function [sys,x0,str,ts]=s_function(t,x,u,flag)
switch flag,
case 0,
[sys,x0,str,ts]=mdlInitializeSizes;
case 1,
sys=mdlDerivatives(t,x,u);
case 3,
sys=mdlOutputs(t,x,u);
case {2, 4, 9 }
sys = [];
otherwise
error(['Unhandled flag = ',num2str(flag)]);
end
function [sys,x0,str,ts]=mdlInitializeSizes
sizes = simsizes;
sizes.NumContStates = 4;
sizes.NumDiscStates = 0;
sizes.NumOutputs = 2;
sizes.NumInputs = 1;
sizes.DirFeedthrough = 0;
sizes.NumSampleTimes = 0;
sys=simsizes(sizes);
x0=[1 0 0 0];
str=[];
ts=[];
function sys=mdlDerivatives(t,x,u)
m=1340;
a=1.54;
b=1.29;
l=a+b;
Ic=2395; %绕质心的转动惯量
rou=sqrt(Ic/m);
k1=40*1000;
k2=44*1000;
M=[m*(b^2+rou^2)/l^2,m*(a*b-rou^2)/l^2;m*(a*b-rou^2)/l^2,m*(a^2+rou^2)/l^2];
K=[k1,0;0,k2];
A=[zeros(2,2),eye(2);-M\K,zeros(2,2)]; %四阶参数4X4矩阵X'=AX
sys=A*x;
function sys=mdlOutputs(t,x,u)
sys(1)=x(1);
sys(2)=x(2);
坐标x1、x2与主坐标x1p、x2p的关系 plot([0;1],[0,0;mode_vec(2),mode_vec(4)]) plot(y(1,1:40),y(2,1:40))
建立系统的动力学方程
建立系统的动力学方程的方法：1、牛顿力学：牛顿第二定律；2、分析力学：拉格朗日方程。 以x、θ为两个变量建立二自由度系统动力学方程；3、影响系数法-张量算法 1、 根据牛顿第二定律，
2、 拉格朗日方程，
d∂L idt∂q
-∂L∂qi
=Qi(i=1,2,.....n) -------------拉格朗日第二类方程
L=T-V，称为拉氏函数，泛函；
势能函数：V=V(q1,q2,.......qn)；
1,q 2,.......q n)； 动能函数：T=T(q
广义坐标qi(i=1,2,.....n)对应的非保守力：Qi(i=1,2,.....n)
d∂L idt∂q
-∂L∂qi
=0(i=1,2,.....n) -------------保守系统的拉氏方程
-------------------利用上诉拉氏方程求解-------------------------------------
T=T=V=
121212
TMq ；V=q(mx
2
2
+Jθ )
12
qKq
T
[k1(x-l1θ)+k2(x+l2θ)]
22
将L=T-V代入拉氏方程可解的：
可见：与牛顿第二定律求得的结果一致。
+Kq=Q;对保守系统Q=0； M q
n
∑(m
j=1
ij
j+kijqj)=Qi(i=1,2,.....n) q
M=(
m0
0J
)
K=(
k1+k2k2l2-k1l1
k2l2-k1l1k1l1+k2l2
2
2
)
q=() θ
j+kijqj)=Qi(i=1,2,.....n)有明确的物理意义：弹性恢复力 +Kq=Q 和∑(mijqM q
j=1n
x
与保守力Q平衡。如张量理论，mij、kij可认为是张量的坐标，kij表-Kq、惯性力-M q
示：使系统仅产生沿qj坐标的单位位移时，沿qi坐标必须施加的外力Qi，或者说Q的i分量在q的j分量上的影响量的投影。
3、影响系数法-张量算法
！！！！！！！下面由张量分量投影计算理论直接求解-----------！！！！！！！！ 设
m11 * + k11 k12 * x
m21 θ + k21 k22 θ
k11:仅当x动1单位时，x向的作用力为k1*1+k2*1=k1+k2;
k12:仅当θ动1单位时，x向的作用力为-k1*(l1*1)+k2*(l2*1)=-k1*l1+k2*l2;
k21:仅当x动1单位时，θ向的作用力为-(k1*1)*l1+(k2*1)*l2=-k1*l1+k2*l2;
k22:仅当θ动1单位时，θ向的作用力为 k1*(l1*1)*l1+k2*(l2*1)*l2=k1*l1^2+k2*l2^2;
m11:仅当x”动1单位时，x”向的作用力为m*1;
k12:仅当θ”动1单位时，x”向的作用力为0; %仅绕质心转动时不影响x”向惯性力 k21:仅当x”动1单位时，θ”向的作用力为0; %仅平动时不影响θ”向惯性力？？？？？？不过质心时该怎么计算？？若旋转中心偏离质心a,则变为
ma
此时，k12: 仅当θ”动1单位时，x”向的作用力为m*(a*1); k21:仅当x”动1单位时，θ”向的作用力为m*1*a;
k22:仅当θ”动1单位时，θ”向的作用力为 (J+m*a^2)*1;
k22:仅当θ”动1单位时，θ”向的作用力为 J*1;
可见与上述结果一致。
对建立的动力学方程更换坐标求偏频
对上述系统建立前后轮纵向位移x1、x2的动力学方程 :x1=x-l1*θ；x2=x+l2*θ。使用matlab的solve（‘x1=x-l1*θ；x2=x+l2*θ’，‘x1’, ‘x2’）可以直接解出: θ=(x1-x2)/(l1+l2);x=(x1*l2+l1*x2)/(l1+l2),代入前面创建的方程组：
消去x1和x2,可得到如下的方程：
ml2+ρl
2
22
2
1+mx
l1l2-ρl
2
2
2+k1x1=0x
m
l1+ρl
2
2
2+mx
l1l2-ρl
2
2
1+k2x2=0x
⎡l22+ρ2⎢m2
l
所以⎢2
⎢ml1l2-ρ
2⎢l⎣
2
2
l1l2-ρ⎤m⎥⎡ 2 1⎤⎡k1xl⨯⎥⎢⎥+⎢22
2⎦⎣0xl+ρ⎥⎣
m12
⎥l⎦
0⎤⎡x1⎤
⎥⨯⎢⎥=0 k2⎦⎣x2⎦
1+η1* 2+ω1*x1=0xx
2
2+η2* 1+ω2*x2=0xx
l1l2-ρ
2
式中，η1=l2+ρ2 联系系数，表示两坐标之间的联系
2η2=
l1l2-ρl1+ρ
2
22
2
2
ω1=
k1l
2
m(l2+ρ)
偏频，表示前后悬挂独立振动时的振动频率，即x1=0时的振动频
率是w2,x2=0时的振动频率是w1,不同于系统的固有频率（2自由度独立时才相等）。 ω2=
22
m(l1+ρ)
k2l
2
ρ=
J
m 汽车绕质心轴的回转半径
在汽车设计中，希望行车时一个悬挂的振动不传到另一个悬挂上，为此，应使车身质量分布和前后轮位置满足：
质量分配系数ε=ll=1，这时η1=η2=0
12ω1=
k1lml2k2lml1
2
ρ
2
=
m1 k2
k1
ω2=
=
m2
2
2
J=mρ
=m1l1+m2l2
ε≠1
对于一般质量分配系数的耦合情况，可以用模态分析法求固有频率及其通解：
x1=A11sin(p1t+ϕ1)+A12sin(p2t+ϕ2)x2=A21sin(p1t+ϕ1)+A22sin(p2t+ϕ2)=β1A11sin(p1t+ϕ1)+β2A12sin(p2t+ϕ2)
⎡x1⎤⎡1即⎢x⎥=⎢β⎣2⎦⎣11⎤⎡A11sin(p1t+ϕ1)⎤⎡1⎥⎢⎥=⎢β2⎦⎣A12sin(p2t+ϕ2)⎦⎣β11⎤⎡xp1⎤
⎥⎥⎢
β2⎦⎣xp2⎦
可见，特征向量阵(模态矩阵)ф组成坐标变换矩阵（由老基到新的主坐标基的坐标变换矩阵），xp=фx为主坐标，主振动的坐标，在该坐标系上，各自由度独立振动（解耦）。
⎡1⎤⎡1⎤
⎢⎥⎢⎥⎣β1⎦⎣β2⎦
T
分别是新的主坐标基的两个基矢量在老基下的投影坐标。
在新的主坐标基下，Mp=фMф为主质量阵（主质量组成的对角阵）；Kp=фKф为主刚度阵（主刚度组成的对角阵）。这种矩阵变换的本质是张量的坐标变换。
+Kpx=0 主坐标方程组为解耦方程组。 Mp x
T
T
利用特征值分解找到系统的主坐标基，通过坐标变换进行解耦、简化计算、然后再变换回去，这是坐标变换的意义所在。
用matlab特征值分解法求平等与转动主模态（振型）
%SH760小轿车空载主要参数 m=1340; a=1.54;
b=1.29;
Ic=2395; %绕质心的转动惯量 k1=40000; k2=44000; M=[m,0;0,Ic];
K=[k1+k2,-(k1*a-k2*b);-(k1*a-k2*b),k1*a^2+k2*b^2];
[eig_vec,eig_val] = eig(inv(M)*K);
[omeg,w_order] = sort(sqrt(diag(eig_val))); %频率 mode_vec = eig_vec(:,w_order); %振型 T=2.*pi./omeg; %周期
mode_vec(:,1)=mode_vec(:,1)./mode_vec(1,1); mode_vec(:,2)=mode_vec(:,2)./mode_vec(1,2); subplot(2,1,1)
plot([1;2],mode_vec(:,1))
title(strcat('w1=',num2str(omeg(1)))); subplot(2,1,2)
plot([1;2],mode_vec(:,2))
title(strcat('w2=',num2str(omeg(2))));
因为对特征值进行了排序，所以w1
求平动与平动主模态（振型）与解析法仿真计算
%SH760小轿车空载主要参数 clear; m=1340; a=1.54; b=1.29;
l=a+b;
Ic=2395; %绕质心的转动惯量 rou=sqrt(Ic/m); k1=40*1000; k2=44*1000;
M=[m*(b^2+rou^2)/l^2,m*(a*b-rou^2)/l^2;m*(a*b-rou^2)/l^2,m*(a^2+rou^2)/l^2]; K=[k1,0;0,k2];
%用matlab特征值分解法求主振型------------------------------------------------
[eig_vec,eig_val] = eig(inv(M)*K);
[omeg] = (sqrt(diag(eig_val))); %频率不用sort排序 mode_vec = eig_vec;%(:,w_order); %振型 T=2.*pi./omeg; %周期
mode_vec(:,1)=mode_vec(:,1)./mode_vec(1,1); mode_vec(:,2)=mode_vec(:,2)./mode_vec(1,2);
w1=sqrt((k1*l^2)/(m*(b^2+rou^2))); w2=sqrt((k2*l^2)/(m*(a^2+rou^2))); w1_pian=sqrt((k1*l)/(m*b)); w2_pian=sqrt((k2*l)/(m*a));
subplot(2,2,1)
plot([1;2],mode_vec(:,1))
title(strcat('w_1=',num2str(omeg(1)),';w_1pian=',num2str(w1_pian))); subplot(2,2,2)
plot([1;2],mode_vec(:,2))
title(strcat('w_2=',num2str(omeg(2)),';w_2pian=',num2str(w2_pian)));
%完全用matlab sym解析法运算求解------------------------------------------------
syms A11 A12 phi1 phi2 t
x1=A11*sin(omeg(1)*t+phi1)+A12*sin(omeg(2)*t+phi2);
x2=mode_vec(2,1)*A11*sin(omeg(1)*t+phi1)+mode_vec(2,2)*A12*sin(omeg(2)*t+phi2); dx1=diff(x1); dx2=diff(x2);
x1_0=subs(x1,'t',0); x2_0=subs(x2,'t',0); dx1_0=subs(dx1,'t',0);
dx2_0=subs(dx2,'t',0);
eq=[sym(strcat(char(x1_0),'=1'));sym(strcat(char(dx1_0),'=0')); sym(strcat(char(x2_0),'=0'));sym(strcat(char(dx2_0),'=0'))];
s=solve_sym(eq);
x1=s.A11(1)*sin(omeg(1)*t+s.phi1(1))+s.A12(1)*sin(omeg(2)*t+s.phi2(1));
x2=mode_vec(2,1)*s.A11(1)*sin(omeg(1)*t+s.phi1(1))+mode_vec(2,2)*s.A12(1)*sin(omeg(2)*t+s.phi2(1)); ti=0:0.02:10;
x1i=subs(x1,'t',ti); x2i=subs(x2,'t',ti);
subplot(2,2,3)
plot(ti',[x1i',x2i'])
%用matlab指数运算求解----------------------------------------------------------
x0=[1;0];xd0=[0;0]; %初始条件 tf=10;dt=0.02; %时间向量
A=[zeros(2,2),eye(2);-M\K,zeros(2,2)]; Y=[x1;x2;x1';x2']
%expm(A)的意义是将坐标先变换到主坐标系，对对角值进行exp运算后再变换到原坐标系，如同张量坐标变换help expm
y0=[x0;xd0]; %四元变量的初始条件
for i=1:round(tf/dt)+1 %设定计算点，作循环计算 tj(i)=dt*(i-1);
y(:,i)=expm(A*tj(i))*y0; %循环计算矩阵指数 end
subplot(2,2,4),plot(tj,[y(1,:)',y(2,:)']),grid
%四阶参数矩阵
Y'=AY-->Y=expm(A*t)*Y0
可见，坐标的选取对固有频率没有影响，但对振型有影响。W1_pianpin和w2_pianpin是前后的偏频（假设质量分配系数为1计算）。
用matlab ode45()直接进行仿真计算
%SH760小轿车空载主要参数 clear; m=1340; a=1.54; b=1.29; l=a+b;
Ic=2395; %绕质心的转动惯量 rou=sqrt(Ic/m); k1=40*1000; k2=44*1000;
M=[m*(b^2+rou^2)/l^2,m*(a*b-rou^2)/l^2;m*(a*b-rou^2)/l^2,m*(a^2+rou^2)/l^2]; K=[k1,0;0,k2];
%用matlab ode45数值解------------------------------------------------
A=[zeros(2,2),eye(2);-M\K,zeros(2,2)]; %四阶参数4X4矩阵X'=AX-->X=expm(A*t)*X0 X=[x1;x2;x1';x2']
syms x1 x2 dx1 dx2
df_sym=A*[x1;x2;dx1;dx2];
df_sym=subs(df_sym,{'x1','x2','dx1','dx2'},{'x(1)','x(2)','x(3)','x(4)'}); n=length(df_sym);
i=1;
ss='['; %先定义好很重要，否则再循环体中定义时，每一循环ss不累加。 while i
ss=strcat(ss,char(df_sym(i)),';'); i=i+1;
end
ss=strcat(ss,char(df_sym(i))); ss=strcat(ss,']'); f=inline(ss,'t','x');
[t,x]=ode45(f,[0 10],[1,0,0,0]);%初始y=0,y'=1 %subplot(2,2,1)
plot(t,[x(:,1),x(:,2)]) %时间状态系列
用s-function进行仿真计算
%sh760.m
function [sys,x0,str,ts]=s_function(t,x,u,flag)
switch flag,
case 0,
[sys,x0,str,ts]=mdlInitializeSizes;
case 1,
sys=mdlDerivatives(t,x,u);
case 3,
sys=mdlOutputs(t,x,u);
case {2, 4, 9 }
sys = [];
otherwise
error(['Unhandled flag = ',num2str(flag)]);
end
function [sys,x0,str,ts]=mdlInitializeSizes
sizes = simsizes;
sizes.NumContStates = 4;
sizes.NumDiscStates = 0;
sizes.NumOutputs = 2;
sizes.NumInputs = 1;
sizes.DirFeedthrough = 0;
sizes.NumSampleTimes = 0;
sys=simsizes(sizes);
x0=[1 0 0 0];
str=[];
ts=[];
function sys=mdlDerivatives(t,x,u)
m=1340;
a=1.54;
b=1.29;
l=a+b;
Ic=2395; %绕质心的转动惯量
rou=sqrt(Ic/m);
k1=40*1000;
k2=44*1000;
M=[m*(b^2+rou^2)/l^2,m*(a*b-rou^2)/l^2;m*(a*b-rou^2)/l^2,m*(a^2+rou^2)/l^2];
K=[k1,0;0,k2];
A=[zeros(2,2),eye(2);-M\K,zeros(2,2)]; %四阶参数4X4矩阵X'=AX
sys=A*x;
function sys=mdlOutputs(t,x,u)
sys(1)=x(1);
sys(2)=x(2);
坐标x1、x2与主坐标x1p、x2p的关系 plot([0;1],[0,0;mode_vec(2),mode_vec(4)]) plot(y(1,1:40),y(2,1:40))