:@o 修改之后我这里能运行啊
b=2.0:0.1:15.0;
options = odeset('RelTol',1e-6,'AbsTol',[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];
[t,X]=ode45('www',ts,x0,options,b(n));
figure(1)
plot(b(n),X(5000:100:10000,1),'.');
xlabel('\fontsize{18}\omega');
ylabel('\fontsize{18}x');grid
hold on
figure(2)
plot(b(n)*30/pi,X(5000:100:10000,2),'.');
xlabel('\fontsize{18}\omega');
ylabel('\fontsize{18}y');grid
hold on
end
function ff=www(t,x,flag,b)
ff=zeros(4,1);
g=9800;
E=0.12;
f=0.12;
a=0.5;
%b=10.8;
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(4,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*t);
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*t)-G;
%ff=[dx(1);dx(2);dx(3);dx(4)] |