下面是程序,谢谢!
麻烦高手们帮忙看看,是什么问题???
function dy=huadian(t,x,flag,w)
m=4;
k=0.25e6;
c=1200;
f=0.2;
q=0.08;
d=0.15;
U=q/d;
w0=sqrt(k/m);
g=9800;
kc=6e7;
e=sqrt(x(1)^2+x(3)^2);
Fx=-d*(1-1/e)*kc*(x(1)-f*x(3))*1/2*(sign(abs(e-1))+sign(e-1));
Fy=-d*(1-1/e)*kc*(x(1)*f+x(3))*1/2*(sign(abs(e-1))+sign(e-1));
dy=[x(2);
Fx/(d*m*w0^2)-c/(m*w0)*x(2)-k*x(1)/(m*w0^2)+w^2*U*cos(w*t);
x(4);
Fy/(d*m*w0^2)-c/(m*w0)*x(4)-k*x(3)/(m*w0^2)+w^2*U*sin(w*t)-g/(d*w0^2)];
clc;
clear all
w=7;
options = odeset('RelTol',1e-5,'AbsTol',[1e-6 1e-6 1e-6 1e-6]);
T=2*pi/w;
ts=[0:T/100:2200*T];
x0=[0.001; 0 ;0.001; 0];
[t,X]=ode45('huadian',ts,x0,options,w);
% % 轴心轨迹图
% figure(1)
% plot(X(2000*100:end,1),X(2000*100:end,3))
% % 庞加莱图
% figure(2)
% plot(X(2000*100:100:end,1),X(2000*100:100:end,3),'.') |