打靶法求弗洛凯特征乘子

飞翼

这是系统随参数u的分岔图
function dbf4
w=1.0;u=0.04;Pm=0.1;Pa=0.18;g=0.1;
T=2*pi/w;h=T/100;m=100;t0=0.0;ep=0.000001;
y0=[0;0];
options=odeset;options.RelTol=1e-6;options.AbsTol=1e-6;
ts=[0:T/100:250*T];
[t,dx]=ode45(@sb,ts,y0,options);%积分250周期作为迭代初值
x(1)=dx(end,1);x(2)=dx(end,2);a=1;b=1;
c(1)=h/2;c(2)=c(1);c(5)=c(1);
c(3)=h;c(4)=h;s(1)=x(1);s(2)=x(2);
v=1;
while v==1
    t=t0;
    x(1)=s(1);x(2)=s(2);
    x1(1)=1.0;x1(2)=0.0;
    x2(1)=0.0;x2(2)=1.0;
    for i=1:4*m  %m表示迭代一周期的次数
        t1=t;ts=t;
        for ii=1:2
            p(ii)=x(ii);w(ii)=x(ii);
        end
        for jj=1:4
            f(1)=p(2);
            if p(1)>=1
            f(2)=-2*u*p(2)-(1+g*cos(t))*(p(1)-1)+Pm+Pa*cos(t);
            elseif p(1)<1&p(1)>-1
            f(2)=-2*u*p(2)++Pm+Pa*cos(t);
            else p(1)<=-1
             f(2)=-2*u*p(2)-(1+g*cos(t))*(p(1)+1)+Pm+Pa*cos(t);
            end
            for ii=1:2
                p(ii)=c(jj).*f(ii)+w(ii);
                x(ii)=c(jj+1).*f(ii)/3.0+x(ii);%在时间间隔h上利用RK算法积分一次
            end
        end
        t=t1;ts=t;
        for ii=1:2
            p(ii)=x1(ii);
            w(ii)=x1(ii);
        end
        for jj=1:4
            f(1)=p(2);
            if x(1)>=1
            f(2)=-(1+g*cos(t))*p(1)-2*u*p(2);
            elseif x(1)<1&x(1)>-1
            f(2)=-2*u*p(2);
            else x(1)<=-1
            f(2)=-(1+g*cos(t))*p(1)-2*u*p(2);
            end
            t=ts+c(jj);
            for ii=1:2
                p(ii)=c(jj).*f(ii)+w(ii);
                x1(ii)=c(jj+1).*f(ii)/3.0+x1(ii);%在时间间隔h上利用RK算法积分一次
            end
        end
        t=t1;ts=t;
        for ii=1:2
            p(ii)=x2(ii);
            w(ii)=x2(ii);
        end
        for jj=1:4
            f(1)=p(2);
            if x(1)>=1
            f(2)=-(1+g*cos(t))*p(1)-2*u*p(2);
            elseif x(1)<1&x(1)>-1
            f(2)=-2*u*p(2);
            else x(1)<=-1
            f(2)=-(1+g*cos(t))*p(1)-2*u*p(2);
            end
            t=ts+c(jj);
            for ii=1:2
                p(ii)=c(jj).*f(ii)+w(ii);
                x2(ii)=c(jj+1).*f(ii)/3.0+x2(ii);%在时间间隔h上利用RK算法积分一次
            end
        end
        b=1;N(a,b)=t;
        b=b+1;N(a,b)=x(1);
        b=b+1;N(a,b)=x(2);
        a=a+1;
        N
    end
    r(1)=s(1)-x(1);
    r(2)=s(2)-x(2);
    dr(1,1)=1.0-x1(1);
    dr(1,2)=-x2(1);
    dr(2,1)=-x1(2);
    dr(2,2)=1.0-x2(2);
    I=[1 0;0 1];
    P=I-dr;%雅可比矩阵
    D=eig(P);%特征乘子
    d=dr(1,1).*dr(2,2)-dr(1,2).*dr(2,1);
    b=-r(2).*dr(1,2)+r(1).*dr(2,2);
    e=-r(1).*dr(2,1)+r(2).*dr(1,1);
    ds(1)=b/d;
    ds(2)=e/d;
    if abs(r(1))<ep&abs(r(2))<ep
        break
    else
        s(1)=s(1)-ds(1);
        s(2)=s(2)-ds(2);
    end
end
D
function dx=sb(t,x)
w=1.0;u=0.04;Pm=0.1;Pa=0.18;g=0.1;
if x(1)>=1;
dx=[x(2);-2*u*x(2)-(1+g*cos(w*t))*(x(1)-1)+Pm+Pa*cos(w*t)];
elseif x(1)<1&x(1)>-1;
   dx=[x(2);-2*u*x(2)+Pm+Pa*cos(w*t)];
else x(1)<=-1;
    dx=[x(2);-2*u*x(2)-(1+g*cos(w*t))*(x(1)+1)+Pm+Pa*cos(w*t)];
end
这是齿轮系统打靶法程序（4周期打靶）
按u从大到小变化
u=0.048时系统开始变为4周期运动
u=0.034时系统开始变为8周期运动
程序结果为：这期间特征乘子都在单位圆内，表示为恒定为4周期运动
与分岔图不符（分岔图明显有分岔行为）；还有就是这期间随u变化，特征乘子还出现共轭复数不止一次，比如：u=0.039时，特征乘子为-0.2197+0.3042i与-0.2197-0.3042i；u=0.038时，特征乘子为：-0.2129+0.2864i与-0.2129-0.2864i；
不知道出现共轭复数正不正常？或是程序有问题，请各位指正！

寒竹冷霜

出现共轭复数是很正常的，楼主我想请教一个问题，困扰我很长时间了，就是你在迭代计算过程中出现过迭代不受敛的情况吗？一般这种情况是由什么原因引起的啊？怎么解决啊，诚求指点啊

飞翼

寒竹冷霜 发表于 2013-11-8 15:26

                               
登录/注册后可看大图

出现共轭复数是很正常的，楼主我想请教一个问题，困扰我很长时间了，就是你在迭代计算过程中出现过迭代不受 ...

确实出现过，但是具体原因我也没搞明白，抱歉，仿真确实很让人头疼

黄忠河

转子动力学里面转子碰摩运动怎么求周期运动的稳定性啊 大神