inoran 发表于 2010-3-25 00:47

ODE解動力三自由度(有程式碼)

我是做一個Mx;''+Kx=P的題目
M=
K=
P=
回圈的方法我參照之前的大大寫的,不過做出來有點問題
所以我就拆開來解六條聯立方程式,內容確認過,應該沒有計算錯誤
的地方!
不過我解出來的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=
yinit=
=ode45('testu',tinterv,yinit)
plot(t,x(:,1))

inoran 发表于 2010-3-25 00:48

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

inoran 发表于 2010-3-25 00:57

我弄成WORD檔了
页: [1]
查看完整版本: ODE解動力三自由度(有程式碼)