单自由度体系刚度随时间变化的振动响应

xinquan

我现在用MATLAB编程计算一个单自由度体系刚度随时间变化的振动响应，激励为余弦波，采用的是WILSON法，程序出来了，结果也能出来，但是我觉得不正确！得到的结果应该是两种谐波的叠加，但是现在出来的结果不是预料的形式，请大家帮忙看下错在何处？谢谢！
程序：
function x=wilsonfa
fid=fopen('wil.txt','w');
m=1e7; c=0; dt=0.01; tao=1.4*dt;
t1=0:0.01:5; k1=6.4e9*ones(1,501); t2=5.01:0.01:20; k2=6.4e9-4e8*(t2-5);
t=[t1 t2]; k=[k1 k2]; f=2*cos(16*pi*t);
df=2*cos(16*pi*(t+tao))-2*cos(16*pi*(t+dt));
xd0=1; xdd0=0;
for i=1:2001
dk(i)=k(i)+(6/tao^2)*m+(3/tao)*c;
ddf(i)=df(i)+m*((6/tao)*xd0+3*xdd0)+c*(3*xd0+(tao/2)*xdd0);
ddx(i)=ddf(i)/dk(i);
ddxd(i)=tao*xdd0+(3/tao)*ddx(i)-3*xd0-(3*tao/2)*xdd0;
ddxdd(i)=(6/tao^2)*ddx(i)-(6/tao)*xd0-3*xdd0;
dxdd(i)=ddxdd(i)/1.4;
dxd(i)=dt*xdd0+(dt/2)*dxdd(i);
dx(i)=dt*xd0+dt^2*0.5*xdd0+(dt^2/6)*dxdd(i);
xd0=xd0+dxd(i);
xdd0=xdd0+dxdd(i);
fprintf(fid,'%10.6f',xdd0);
end
fid2=fopen('wil.txt','rt');
x=fscanf(fid2,'%f',[1,2001]);
plot(x)

[ 本帖最后由 ChaChing 于 2009-7-12 12:16 编辑 ]

lxq

那些数据文件没有没法运行呀？<BR>帮你顶一下~~

cdwxg

瓜啊你<BR>人家说了程序出来了，结果也出来了<BR>但不是预料的那样，所以认为程序有错...不过我也不会谐波的东西<BR>程序嘛哪里错了还可能行<BR>要问原理类的东西，这个方面我不行。