马上注册,结交更多好友,享用更多功能,让你轻松玩转社区。
您需要 登录 才可以下载或查看,没有账号?我要加入
x
计算出的x0不为0,为什么画出的图像值为0哪 ?想了好几天不知道为什么,哪位大侠能帮忙看下 不胜感激!!
%Wilson - theta 法
function wilson(m ,c ,k ,x0 ,v0 ,time ,delt)
m=[60000 0;0 50000]; k=[80000 -30000;-30000 30000]; c=[3.5 -1.5;-1.5 1.5];
x0=[0;0]; v0=[0;0]; time=100; delt=0.1;
n = time/ delt ; disp = zeros(3 ,n) ; %设定存储位移矩阵
m_inv = inv(m + delt*c/2 +(delt^2) * k/6);
[mod ,fre ] = eig(inv(m)*k) ; diag(sqrt(fre) ) ; %固有频率
i = 1 ; beta = 1.4 ;
for t = 0:delt:time
f = [2.0 * sin(25 * t);-3.0 *cos(2.2 *t)]; %简谐作用力向量,需要根据实际情况确定输入
if t == 0, xdd0=inv(m)*(f-k*x0-c*v0); %初始加速度
else
xdd=m_inv*(f-c*(v0 +(delt/2)*xdd0)-k*(x0 + delt*v0 + (1/2 - beta) *delt^2 *xdd0) ) ;
x=m_inv*(m*(x0 + delt*v0 +delt^2/ 3*xdd0) + c*(delt/ 2*x0 + delt^2/3* v0 + delt^3/ 12*xdd0) + delt^2/ 6*f) ; %计算位移
xd=v0+delt*(xdd0 + xdd)/2 ; %计算速度
xdd0 = xdd ; v0 = xd ; x0 = x ;
disp( ': ',i) = x0; %显示质点位移
i = i + 1 ;
end
end
t = 1:n ;
figure ('numbertitle','off','name','第一层框架位移','pos',[450 180 400 420 ]) ;
plot (t,disp ( 1,:)) ; grid ;
[ 本帖最后由 ChaChing 于 2009-3-5 23:33 编辑 ] |