请教：模态叠加法和Runge-Kutta方法解动力学方程的区别

tianlei2005

一个常见的三自由度质量－弹簧系统，其动力学方程为：
[M]{x''}+[K]{x}={F}
质量、刚度和激励矩阵分别为：
M=diag([1;1;1]);k=[3 -1 0;-1 2 -1;0 -1 3];F={sin(3*t);0;0};
我分别用模态叠加法和Runge-Kutta算法求解，但是两种解法得到的结果却不相同，请问这是什么原因，何种方法才是正确的。
这是我用模态叠加法的代码：
%利用模态叠加法分析扭转振动的实例
%定义符号,a、b为阻尼的比例系数
syms t;
%齿轮系统参数
m1=1;m2=1;m3=1;k1=2;k2=2;k3=1;k4=2;
M=diag([m1;m2;m3]);
% K=[k1+k2 -k2 0;
%     -k2 k2+k3 -k3;
%     0 -k3 k3+k4];
K=[3 -1 0;-1 2 -1;0 -1 3];
F1=1;F2=0;F3=0;w=3;
w0(1)=w;w0(2)=0*w;w0(3)=0*w;
f1=F1*sin(w0(1)*t);f2=F2*sin(w0(2)*t);f3=F3*sin(w0(3)*t);
%求振动系统的正则振型矩阵和固有频率
[N,V]=eig(K,M);
ModeValue=sqrt(diag(V));
w1=ModeValue(1);
w2=ModeValue(2);
w3=ModeValue(3);
%定义振动系统的初始条件
x0=[0 0 0]';
X0=N'*M*x0;
    a1=X0(1);
    a2=X0(2);
    a3=X0(3);
%正则变换，求出正则坐标下的激励力
Q=[f1;f2;f3];
Q_Fai=N'*[f1;f2;f3];
Qn=N'*diag([F1;F2;F3]);
% Qn1=N'*[f1;0;0];
% Qn2=N'*[0;f2;0];
% Qn3=N'*[0;0;f3];
Q_Fai1=N'*[f1;0;0];
Q_Fai2=N'*[0;f2;0];
Q_Fai3=N'*[0;0;f3];
%求出在正则坐标下的响应
n=length(M);
%求系统的无阻尼受迫振动解
for i=1:3
    b1(i)=1/(ModeValue(i)^2-w^2);
    b2(i)=1/(ModeValue(i)^2-(3*w)^2);
    b3(i)=1/(ModeValue(i)^2-(8*w)^2);
end
Xn1=b1'.*Q_Fai1;
Xn2=b2'.*Q_Fai2;
Xn3=b3'.*Q_Fai3;
x1=N*Xn1;
x2=N*Xn2;
x3=N*Xn3;
x=x1+x2+x3;
% x=C*Xn;
% Qn=N'*[F1;F2;F3];
x=-0.088.*[1;2;1].*sin(3*t)-2.63.*[-1.732;0;1.732].*sin(3*t)+0.21.*[1.414;-1.414;1.414].*sin(3*t)
%定义仿真时间和采样点数
t=0:0.01:100;
%对结果进行fft变换
u1=eval(x(1));
u2=eval(x(2));
u3=eval(x(3));（省去后面的画图和FFt变换部分），


下面是调用ode45函数的代码
%test4.m
function f=test4(t,y);
m1=1;m2=1;m3=2;k=1;
U=[0 1 0 0 0 0;
    -3 0 1 0 0 0;
    0 0 0 1 0 0;
    1 0 -2 0 1 0;
    0 0 0 0 0 1;
    0 0 1 0 -3 0];
f=U*y+[0 sin(3*t) 0 0 0 0]';

%test4Result.m
[t,y]=ode45('test4',[0:0.01:100],[0 0 0 0 0 0]');
u1=y(:,1);
u2=y(:,3);
u3=y(:,5);
(后面省去画图和fft变换部分)

21172485

他们是等效的,模态叠加一个是解析方法,可以看清解的结构,
RK是数值方法,多用来分析系统的动力学行为等

wwg.no.19

感觉一下这段好像有问题，

for i=1:3
    b1(i)=1/(ModeValue(i)^2-w^2);可疑
    b2(i)=1/(ModeValue(i)^2-(3*w)^2);
    b3(i)=1/(ModeValue(i)^2-(8*w)^2);
end
Xn1=b1'.*Q_Fai1;可疑
Xn2=b2'.*Q_Fai2;
Xn3=b3'.*Q_Fai3;

解偶后求解单个方程时公式好像用错了，应该时Bd=f0/sqrt(k-m*w^2),而且相位也没有确定啊，鄙人刚学不久，也许指证有错，请指教

欧阳中华

.
  模态叠加法计算时，没有初始条件问题，或不受初始条件影响；而时域分析时，需考虑加载时初始条件，也就是瞬态和稳态的差别，如果模态叠加法与时域的RK比，就要取稳态部分，那么判断稳态的时间段.. ..

vib

我记得振动基础里面有模态叠加法，不知道是不是一个模态叠加法，数值分析里面学的Rugge-Kutta是用来解一个有初值的一维微分方程，请推荐本这方面的书吧，谢谢！

vib

如果没有阻尼的话，象这样的方程，应该可以得到精确解吧，没必要求数值解吧

fineshang

请教各位高手，12自由度的二阶微分方程组的频域特性曲线借助matlab，用四阶龙格库塔求解出了时域解，那接下来的频域特性曲线怎么求啊？？
知道的就麻烦告一下吧，非常感谢！！！很着急呀！！！！！！

pengweicai

频域特性比较类似。
首先要把动力学方程转换到频域，后面的就基本类似了。

fineshang

方程形如：[m]{X''}+[C]{X'}+[K(t)]{X}={P}
其中[m]=diag{m1,m2,m3,m4,m5,m6,m7,m8,m9,m10,m11,m12}
    {X}=（x1，x2，x3，x4，x5，x6，x7，x8，x9，x0，x11，x12）
    [C]为12×12阶的常数矩阵
    [K]为12×12阶的时变矩阵，其中每个元素均是时间t的函数，Kij=Kij(t)
    {P}={p1(t),p2(t),p3(t),p4(t),p5(t),p6(t),p7(t),p8(t),p9(t),p10(t),p11(t),p12(t)}’
想用数值解法（四阶龙格库塔法），借助matlab进行求解，最后输出各自变量的频域解！可以解出频域解么？？？？？！非常感谢！知道的或做过这方面研究的朋友，高手，还请指教！！！

fineshang

也就是说每个都进行拉普拉斯变化么？请这位朋友明示，怎么个转换到频域法？谢谢了！！！

fineshang

请问，这位朋友，怎么把动力学方程转换到频域呢？？？
请指教啊，非常感谢！！！！！！！