|
马上注册,结交更多好友,享用更多功能,让你轻松玩转社区。
您需要 登录 才可以下载或查看,没有账号?我要加入
x
function dydt=rk_beamone(t,y,flag,f,c1,c2,c3,k1,k2,cf,m1,m2,A)
w=2*pi*f;
g1=A*sin(w*t);
g2 =w*A*cos(w*t);
t
dydt=[y(2);-(k1*(y(1)-g1)+c1*(y(2)-g2)+k2*(y(1)-y(3))^3+c2*(y(2)-y(4))+cf*sign(y(2)-y(4))+c3*(y(2)-y(4))^3)/m1;
y(4);-(k2*(y(3)-y(1))^3+c2*(y(4)-y(2))+cf*sign(y(4)-y(2))+c3*(y(4)-y(2))^3)/m2];
clc
clear
% h_opt=odeset('RelTol',1e-8,'AbsTol',1e-8,'MaxStep',1e-7);
h_opt=odeset;
m1=60;
k1=2.1346e6;
c1=10;
m2=6;
c2=600;
k2=5e18;
c3=800;
cf=0.48;
gap=1;
A=0.001;
tspan=[0:0.01:2];
% f=20;
h=1;
fmax=50;
n=fmax/h+1;
ymax=zeros(n,2);
i=0;
for f=0:h:fmax
i=i+1;
y1=[0;0.001;0;0];
[t,y0]=ode45('rk_beamone',tspan,y1,h_opt,f,c1,c2,c3,k1,k2,cf,m1,m2,A);
end
% plot(t,y0(:,1),'k-.','LineWidth',2)
ymax(i,1)=max(y0(:,1));
f1=0:h:fmax;
plot(f1,ymax)
|
|