打靶法求Floquet特征乘子

飞翼

function dbf3
m=20;h=0.3142;t0=0.0;ep=0.00000001;
x(1)=0.06594;x(2)=0.01;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);
x1(1)=1.0;x1(2)=0.0;
x2(1)=0.0;x2(2)=1.0;
v=1;
while v==1
    t=t0;
    x(1)=s(1);x(2)=s(2);
    iii=iii+1;
    for i=1: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);
            f(2)=-0.2*p(2)-4.0*p(1)-p(1).^3+0.03*cos(t);
            t=ts+c(jj);
            for ii=1:2
                p(ii)=c(jj).*f(ii)+w(ii);
                x(ii)=c(jj+1).*f(ii)/3.0+x(ii);
            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);
            f(2)=-0.2*p(2)-(4.0+3.0*x(1).^2).*p(1);
            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);
            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);
            f(2)=-0.2*p(2)-(4.0+3.0*x(1).^2).*p(1);
            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);
            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;
    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(ds(1))<ep&abs(ds(2))<ep
        break
    else
        s(1)=s(1)-ds(1);
        s(2)=s(2)-ds(2);
    end
end
N
P
D
这是打靶法求达芬系统的1T周期解的程序，我同时加了点同求系统的雅可比矩阵P和其特征值D，得到的结果是
P =
  1.0e-005 *
   0.132237285277448   0.049603766980269
  -0.197319424225645   0.141870719549964

D =
  1.0e-005 *
  0.137054002413706 + 0.098815919648602i
  0.137054002413706 - 0.098815919648602i
发现值都非常小，不知道合不合理，也不知道这么求对不对
另外有点疑惑的是蓝色标注行，dr(1,2)=x2(1);dr(2,1)=x1(2);是不是x2（1），x1（2）前应该加负号
b=-r(2).*dr(1,2)+r(1).*dr(2,2);e=-r(1).*dr(2,1)+r(2).*dr(1,1);是不是应该改为b=r(2).*dr(1,2)-r(1).*dr(2,2);e=r(1).*dr(2,1)-r(2).*dr(1,1);
请各位帮忙看下
   

寒竹冷霜

感觉dr(1,2)=x2(1);dr(2,1)=x1(2);这里应该不仅要加负号还要换个位置
dr(1,2)=-x1(2);dr(2,1)=-x2(1)
后面的b和d也不对劲
不知道对不对啊

飞翼

寒竹冷霜 发表于 2013-9-11 20:49

                               
登录/注册后可看大图

感觉dr(1,2)=x2(1);dr(2,1)=x1(2);这里应该不仅要加负号还要换个位置
dr(1,2)=-x1(2);dr(2,1)=-x2(1)
后面 ...

应该不是吧，dr（1,2）表示r（1）对s（2）求导
dr（2,1）表示r（2）对s（1）求导

寒竹冷霜

飞翼 发表于 2013-9-13 08:38

                               
登录/注册后可看大图

应该不是吧，dr（1,2）表示r（1）对s（2）求导
dr（2,1）表示r（2）对s（1）求导

你这是对的，但是r（1）对s（2）求导不应该是x1(2)吗？
不过这都不影响，这是个对称阵，问题在于加上负号后，b和e也加负号
怎么迭代就进行不下去啊

飞翼

寒竹冷霜 发表于 2013-9-13 10:01

                               
登录/注册后可看大图

你这是对的，但是r（1）对s（2）求导不应该是x1(2)吗？
不过这都不影响，这是个对称阵，问题在于加上负号 ...

r（1）对s（2）求导是x2（1）……

寒竹冷霜

飞翼 发表于 2013-9-25 15:56

                               
登录/注册后可看大图

r（1）对s（2）求导是x2（1）……

不知道你在利用打靶法的同时有没有结合continuation technique或者叫path-following technique来算算分叉图之类的

飞翼

寒竹冷霜 发表于 2013-9-30 10:25

                               
登录/注册后可看大图

不知道你在利用打靶法的同时有没有结合continuation technique或者叫path-following technique来算算分叉 ...

这个倒是没有，帮不了你

gclx

可以联系我，18678378320

补充内容 (2013-11-4 15:59):
如果你的雅克比矩阵代表的不是周期初和周期末值之间的映射关系，那么floquet乘子就不是这个雅克比矩阵的特征值。比如某齿轮振动系统方程如下：
     x1=x2;
     x2=-2*jieta*x2-k*x1-p(wt);
他的雅克比矩阵是反应的连续时间下的前一时刻到后一时刻的映射，并不是周期初到周期末的映射，所以这个雅克比矩阵的特征值就不是对应的floquet乘子。


补充内容 (2013-11-4 16:00):
但是对于大多数非线性系统是不能直接写出周期初到周期末的映射的，而需要用数值计算方法将一个周期T分成很多小段，然后一段一段的数值积分得到周期末的值。那么对于这种系统求解floquet乘子的比较笨的方法可以这样：
将【0 T】分成等间隔时间，t1，t2，.....T。从t1时刻计算到t2时刻对以一个矩阵M1，t2到t3时刻对应一个矩阵M2，一次往后，然后求解M=M1*M2*.....Mn，就是状态转移矩阵。然后求解这个状态转移矩阵的特征值就是floquet乘子。步长越小，floquet乘子越精确。数值计算选用方法精度越高，floquet乘子计算越准确。


补充内容 (2013-11-4 16:00):
这时可以把系统一个周期分成t1，t2.....tn,求出每一时刻的矩阵Mi，这个可以用matlab的符号矩阵来处理，最终求出M。
y(t1)=M0*y(t0);
y(t2)=M1*M0*y(t0);
.............
y(tm)=G*y(t0)
那么G就是以上M的累乘积，这个G除了与系统参数有关外，还是初始时间有关系，这个初始时间我想可以直接代入0。那么初值到底取什么y0，G是不会涉及的。

黄忠河

我也在做这个 看不懂啊

黄忠河

谁可以告诉我