一個常見的三自由度質量－彈簧系統，其動力學方程為：
[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變換部分)