Lyapunov exponent for fractional order differential equation

Question

ABHIJIT MAJUMDER 2020-5-4

1
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/523006-lyapunov-exponent-for-fractional-order-differential-equation

编辑： Walter Roberson 2022-7-26

Hi, I have three dimensional fde model of which I want to compute Lyapunov exponent with respect to parameter. I have taken this code from the paper entitled 'Matlab code for Lyapunov exponents of fractional order systems' but it's not working, Please help.

function [t,LE]=FO_Lyapunov_q(ne,ext_fcn,t_start,h_norm,t_end,x_start,h,q,p);
x=zeros(ne*(ne+1),1);
x0=x;
c=zeros(ne,1);
gsc=c; zn=c;
n_it = round((t_end-t_start)/h_norm);
x(1:ne)=x_start;
i=1;
while i<=ne
    x((ne+1)*i)=1.0;
    i=i+1;
end
t=t_start;
it=1;
while it<=n_it
    [T,Y] = fde12(q,ext_fcn,t,t+h_norm,x,h);
    t=t+h_norm;
    Y=transpose(Y);
    x=Y(size(Y,1),:);
    i=1;
    while i<=ne
        j=1;
        while j<=ne;
            x0(ne*i+j)=x(ne*j+i);
            j=j+1;
        end;
        i=i+1;
    end;
    
    zn(1)=0.0;
    j=1;
    while j<=ne
        zn(1)=zn(1)+x0(ne*j+1)*x0(ne*j+1);
        j=j+1;
    end;
    zn(1)=sqrt(zn(1));
    j=1;
    while j<=ne
        x0(ne*j+1)=x0(ne*j+1)/zn(1);
        j=j+1;
    end
    j=2;
    while j<=ne
        k=1;
        while k<=j-1
            gsc(k)=0.0;
            l=1;
            while l<=ne;
                gsc(k)=gsc(k)+x0(ne*l+j)*x0(ne*l+k);
                l=l+1;
            end
            k=k+1;
        end
        k=1;
        while k<=ne
            l=1;
            while l<=j-1
                x0(ne*k+j)=x0(ne*k+j)-gsc(l)*x0
                (ne*k+l);
                l=l+1;
            end
            k=k+1;
        end;
        zn(j)=0.0;
        k=1;
        while k<=ne
            zn(j)=zn(j)+x0(ne*k+j)*x0(ne*k+j);
            k=k+1;
        end
        zn(j)=sqrt(zn(j));
        k=1;
        while k<=ne
            x0(ne*k+j)=x0(ne*k+j)/zn(j);
            k=k+1;
        end
        j=j+1;
    end
    % update running vector magnitudes
    k=1;
    while k<=ne;
        c(k)=c(k)+log(zn(k));
        k=k+1;
    end;
    % normalize exponent
    k=1;
    while k<=ne
        LE(k)=c(k)/(t-t_start);
        k=k+1;
    end
    i=1;
    while i<=ne
        j=1;
        while j<=ne;
            x(ne*j+i)=x0(ne*i+j);
            j=j+1;
        end
        i=i+1;
    end;
    x=transpose(x);
    it=it+1;
end
%----------------------------------------------------------------------------------------------------------------------------%
    function f=LE_RF(t,x,p)
        f=zeros(size(x));
        X= [x(4), x(7), x(10);
            x(5), x(8), x(11);
            x(6), x(9), x(12)];
        %RF equations
        f(1)=x(2)*(x(3)-1+x(1)*x(1))+0.1*x(1);
        f(2)=x(1)*(3*x(3)+1-x(1)*x(1))+0.1*x(2);
        f(3)=-2*x(3)*(p+x(1)*x(2));
        %Jacobian matrix
        J=[2*x(1)*x(2)+0.1, x(1)*x(1)+x(3)-1, x(2);
            -3*x(1)*x(1)+3*x(3)+1,0.1,3*x(1);
            -2*x(2)*x(3),-2*x(1)*x(3),-2*(x(1)*x(2)+p)];
        %Righthand side of variational equations
        f(4:12)=J*X; % To be modified if ne>3
        %----------------------------------------------------------------------------------------------------------------------------%
        function run_Lyapunov_p(ne,ext_fcn,t_start,h_norm,t_end,x_start,h,q,p_min,p_max,n);
            hold on;
            p_step=(p_max-p_min)/n
            p=p_min;
            while p<=p_max
                [t,LE]=FO_Lyapunov_q(ne,ext_fcn,t_start,h_norm,t_end,x_start,h,q,p);
                p=p+p_step
                plot(p,LE);
            end
        end
        %----------------------------------------------------------------------------------------------------------------------------%
        run_FO_Lyapunov_q(3,@LE_RF,0,0.05,150,[0.1;0.1;0.1],0.002,0.9,1,800)