动力板砖 发表于 2010-5-26 11:56

求助高手帮忙看下程序

我做了个自治系统的分岔图,程序是参考论坛中的程序修改的 但是一有问题 请高手帮忙看一下

%%函数主程序
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);
    % 舍弃前面迭带的结果,用后面的结果画图
=ode23('fric',,);
    =ode23('fric',,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=;
            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=;   
            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
    =ode23('fric',,Y(length(Y),:));

happy 发表于 2010-5-26 17:37

ode45 改为 ode15s 看看
页: [1]
查看完整版本: 求助高手帮忙看下程序