|
马上注册,结交更多好友,享用更多功能,让你轻松玩转社区。
您需要 登录 才可以下载或查看,没有账号?我要加入
x
有没有做碰撞问题的,我下面作了一个单自由度的碰撞的matlab程序,但是运行不出结果,不知道问题出在那里,请指教。谢谢!
%方程
%d2x/dt2+k*dx/dt+x=cos(w*t)
%边界条件 x(t+)=x(t-),dx/dt(t+)=-r* dx/dt(t-)
clear;
global w;
a=0.01;
pi=3.1415;
r=0.6;
w=2.5;
fid= fopen('xjl.txt','w');
%fid1= fopen('xjl1.txt','w');
x=1.6;
v=0;
t=0;
t1=0.05
n=1;
while j<500
t=t+j*t1
A=(1-w*w)/((1-w*w)^2+(2*w*a)^2);
B=2*w*a/((1-w*w)^2+(2*w*a)^2);
w1=sqrt(1-a^2);
c=2*n*pi/w;
d=exp(-c*a);
e=-r*d;
w2=w1*c;
f=d*sin(w2)/(d-1);
g=w1+e*a*sin(w2)-e*w1*cos(w2);
c1=f*(e*g*w-g*w)/(g-a*f*cos(w2)-e*f*w1*sin(w2));
c2=(e*g*w-g*w)/(g-a*f*cos(w2)-e*f*w1*sin(w2));
x=exp(-a*t)*(c1*cos(w1*t)+c2*sin(w1*t))+A*sin(w*t)+B*cos(w*t);
v=exp(-a*t)*((c2*w1-a*c1)*cos(w1*t)-(a*c2+c1*w1)*sin(w1*t))+A*w*cos(w*t)-B*w*sin(w*t);
if x<=0
x=x;
v=-r*v;
end
if j>200
fprintf(fid,'%12.8f %12.8f %12.8f \n',t,x,v);
end
end
fclose(fid);
load xjl.txt;
figure
plot(xjl(:,2),xjl(:,3)); |
|