|
马上注册,结交更多好友,享用更多功能,让你轻松玩转社区。
您需要 登录 才可以下载或查看,没有账号?我要加入
x
被调用的函数
function ss=Rugkut3(t,y,flag,A,B,u )
ss=A*y+B*u;
运行的程序代码
clf
clear
clc;
kaxi=0.01; mu=0.02;
kaxit=0.002;
h=0.1;
NT=50;
t=0:h:NT;
N=length(t);
ft=sqrt(1+0.5*mu)/(1+mu);
dalt=0.6;
fe=1.0;
t0=0;
e=0.6;
zz0(:,1)=[0;0;0;0];
for I=1:N
Xg(I)=cos(fe*t(I));
if abs(zz0(2,1)-zz0(1,1))-dalt>=1*0.0000001
A=[0 0 1 0;0 0 0 1;0 mu*ft.^2 -2*kaxi 2*mu*kaxit*ft*e;0 -(1+mu)*ft.^2 2*kaxi
-(1+mu)*2*kaxit*ft*e];
B=[0;0;1;0];
u=Xg(I);
k1=Rugkut3(t0,zz0,A,B,Xg(I));
k2=Rugkut3(t0+h/2,zz0+(h/2)*k1,A,B,Xg(I));
k3=Rugkut3(t0+h/2,zz0+(h/2)*k2,A,B,Xg(I));
k4=Rugkut3(t0+h,zz0+h*k3,A,B,Xg(I));
zz1=zz0+(1/6)*h*(k1+2*k2+2*k3+k4);
else
zz1(1,1)=zz0(1,1);
zz1(2,1)=zz0(2,1);
zz1(3,1)=(e+mu)/(e*(1+mu))*zz0(3,1)+mu*(e+1)/(e*(1+mu))*zz0(4,1);
zz1(4,1)=(1+e)/(e*(1+mu))*zz0(3,1)-(mu*e-1)/(e*(1+mu))*zz0(4,1);
end
Z12(:,I)=zz1(:,1);
t0=t0+h;
zz0=zz1;
end
figure(1)
plot(t,(Z12(1,:)));
title({'系统位移-时间关系图'})
xlabel('t')
ylabel(' x_1')
hold on;
figure(3),
plot(Z12(1,:),Z12(3,:))
title({'系统位移-速度关系图'})
xlabel(' x_1')
ylabel('dx_1/dt')
hold on;
figure(2),
plot(t,Z12(1,:)+dalt,'g:',t,(Z12(2,:)+Z12(1,:)),'r-',t,Z12(1,:)-dalt,'g: ');
title({'碰撞效果图'})
xlabel('t')
ylabel('x_1-\delta x_2 x_1+\delta')
hold on
我刚学matlab,有很多问题不懂,望各位同学能指点指点,我急需要解决这个问题,谢谢!!!!!!!!!! |
|