|
楼主 |
发表于 2007-8-20 17:14
|
显示全部楼层
回复 #86 咕噜噜 的帖子
b=2.0:0.1:20;
options = odeset('RelTol',1e-6,'AbsTol',[1e-6 1e-6 1e-6 1e-6 1e-6]);
for n=1:length(b);
w=2.0;
T=2*pi/w;
ts=[0:T/100:100*T];
x0=[0.001 0 0.001 0 0];
[t,X]=ode45('www',ts,x0,options,b(n));
figure(1)
plot(b(n),X(5000:100:end,1),'k.');
xlabel('\fontsize{18}\omega');
ylabel('\fontsize{18}x');
hold on
end
grid
function ff=www(t,x,flag,b)
ff=zeros(5,1);
g=9800;
E=0.12;
f=0.12;
a=0.5;
%b=15.0;
q=0.16;
f0=25;
u=0.136/q;
w=2.0;
e=sqrt(x(1)^2+x(2)^2);
G=g/((2*pi*f0)^2*q);
ff=zeros(5,1);
dx(1)=x(3);
dx(2)=x(4);
dx(3)=-2*E*x(3)-x(1)-a*(x(1)^2+x(2)^2)*x(1)-b*(1-1/e)*(x(1)-f*x(2))+u*w^2*cos(w*x(5));
dx(4)=-2*E*x(4)-x(2)-a*(x(1)^2+x(2)^2)*x(2)-b*(1-1/e)*(f*x(1)+x(2))+u*w^2*sin(w*x(5))-G;
dx(5)=1;
ff=[dx(1);dx(2);dx(3);dx(4);dx(5)]; |
|