马上注册,结交更多好友,享用更多功能,让你轻松玩转社区。
您需要 登录 才可以下载或查看,没有账号?我要加入
x
我做了个自治系统的分岔图,程序是参考论坛中的程序修改的 但是一有问题 请高手帮忙看一下
%%函数主程序
function dy=fric(t,y)
global v;
%-----参数取值-------
m=0.01;c1=10;c2=1;k1=100000;k2=200000;theta=0.1;sigma0=100000;sigma2=0.01;
miuk=1.0;mius=1.5;vs=1;
%-----系统参数的计算------
cc1=(c1*theta^2+c2)/(m+m*theta^2);cc2=(c2*theta-c1*theta)/(m+m*theta^2);
kk1=(k1*theta^2+k2)/(m+m*theta^2);kk2=(k2*theta-k1*theta)/(m+m*theta^2);
gv=miuk+(mius-miuk)*exp(-(v/vs)^2);
%gv=1.5-167.5042*v;
%Dgv=(2*v/vs^2)*exp(-(v/vs)^2);
Dgv=(v/vs^2)*exp(-(v/vs)^2);
%-----系统方程-------
dy=zeros(3,1);
dy(1)=y(2);
dy(2)=-(kk1+kk2*(gv+sigma2*v))*y(1)-(cc1+cc2*(gv+sigma2*v))*y(2)+kk2*sigma2*y(1)*y(2)-kk2*sigma0*y(1)*y(3)-cc2*sigma0*y(2)*y(3)+cc2*sigma2*y(2)^2;
dy(3)=-v*Dgv/gv*dy(2)-sigma0*v/gv*dy(3);
%%%分岔图程序
function fric_bifur_v
global v
Z=[];
for v=linspace(0.01,10,50);
% 舍弃前面迭带的结果,用后面的结果画图
[T,Y]=ode23('fric',[0,1],[0.0001;0.5;0.00000001]);
[T,Y]=ode23('fric',[0,5],Y(length(Y),:));
Y(:,1)=Y(:,2)-Y(:,1);
% 对计算结果进行判断,如果点满足x=y,则取点
for k=2:length(Y)
f=k-1;
if Y(k,1)<0
if Y(f,1)>0
y=Y(k,2)-Y(k,1)*(Y(f,2)-Y(k,2))/(Y(f,1)-Y(k,1));
Z=[Z v+abs(y)*i];
end
else
if Y(f,1)<0
y=Y(k,2)-Y(k,1)*(Y(f,2)-Y(k,2))/(Y(f,1)-Y(k,1));
Z=[Z v+abs(y)*i];
end
end
end
end
plot(Z,'.','markersize',1)
title('映射分岔图')
xlabel('v'),ylabel('|y| where x=y')
下面是错误提示
Warning: Failure at t=4.354267e-018. Unable to meet integration tolerances without reducing the step size below the smallest value allowed (1.232595e-032) at time t.
> In ode23 at 369
In fric_bifur_v at 9
Error in ==> ode23 at 495
idx = oldnout+1:nout;
Error in ==> fric_bifur_v at 9
[T,Y]=ode23('fric',[0,5],Y(length(Y),:)); |