声振论坛

 找回密码
 我要加入

QQ登录

只需一步,快速开始

查看: 3692|回复: 3

[综合讨论] 打靶法求弗洛凯特征乘子

[复制链接]
发表于 2013-8-25 12:17 | 显示全部楼层 |阅读模式

马上注册,结交更多好友,享用更多功能,让你轻松玩转社区。

您需要 登录 才可以下载或查看,没有账号?我要加入

x
ww.jpg
这是系统随参数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;
不知道出现共轭复数正不正常?或是程序有问题,请各位指正!


QQ图片20130825120145.jpg
回复
分享到:

使用道具 举报

发表于 2013-11-8 15:26 | 显示全部楼层
出现共轭复数是很正常的,楼主我想请教一个问题,困扰我很长时间了,就是你在迭代计算过程中出现过迭代不受敛的情况吗?一般这种情况是由什么原因引起的啊?怎么解决啊,诚求指点啊
 楼主| 发表于 2013-11-12 10:38 | 显示全部楼层

确实出现过,但是具体原因我也没搞明白,抱歉,仿真确实很让人头疼
发表于 2019-4-15 16:29 | 显示全部楼层
转子动力学里面转子碰摩运动怎么求周期运动的稳定性啊 大神
您需要登录后才可以回帖 登录 | 我要加入

本版积分规则

QQ|小黑屋|Archiver|手机版|联系我们|声振论坛

GMT+8, 2024-12-23 06:13 , Processed in 0.067722 second(s), 22 queries , Gzip On.

Powered by Discuz! X3.4

Copyright © 2001-2021, Tencent Cloud.

快速回复 返回顶部 返回列表