ODE解動力三自由度（有程式碼）

inoran

我是做一個Mx;''+Kx=P的題目
M=[m1,0,0, ; 0 ,m2 0 ; 0, 0 ,m3]
K=[k1+k2,-k2,0; -k2, k2+k3,-k3 ; 0, -k3, k3+k4]
P=[Po*sin(w*t);0;0]
回圈的方法我參照之前的大大寫的,不過做出來有點問題
所以我就拆開來解六條聯立方程式,內容確認過,應該沒有計算錯誤
的地方！
不過我解出來的X1（第一個廣義座標的位移）的圖跟講義上的差了許多
我也不知道我錯在那邊！！  謝謝了
副程式：
function dx=testu(t,x)
global w a11 a12 a21 a22 a23 a32 a33 m1 m2 m3 k1 k2 k3 k4
a11=-m1*(k1+k2);
a12=-m1*(-k2);
a21=-m2*(-k2);
a22=-m2*(k2+k3);
a23=-m2*(-k3);
a32=-m3*(-k3);
a33=-m3*(k3+k4);
dx=zeros(6,1);
dx(1)=x(2);
dx(2)=a11*x(1)+a12*x(3)+m1*sin(w*t)
dx(3)=x(4);
dx(4)=a21*x(1)+a22*x(3)+a23*x(5);
dx(5)=x(6);
dx(6)=a32*x(3)+a33*x(5);

主程式：
lear all
global w a11 a12 a21 a22 a23 a32 a33 m1 m2 m3 k1 k2 k3 k4
w=3
m1=1 ;
m2=1 ;
m3=1;
k1=2;
k2=2;
k3=1;
k4=2;
a11=-m1*(k1+k2);
a12=-m1*(-k2);
a21=-m2*(-k2);
a22=-m2*(k2+k3);
a23=-m2*(-k3);
a32=-m3*(-k3);
a33=-m3*(k3+k4);

tinterv=[0 50]
yinit=[0 0 0 0 0 0 ]
[t,x]=ode45('testu',tinterv,yinit)
plot(t,x(:,1))

inoran

正解的圖片他不給我上傳＠＠

inoran

我弄成WORD檔了