关于duffing系统求Lyapunov指数的程序，请各位高手指点

脚大老田

function dX = DF(t,X,flag,alf)global alf;a=1;b=0.25;c=1;omega=1;x=X(1);y=X(2);z=X(3);Y = [X(4), X(7), X(10);    X(5), X(8), X(11);    X(6), X(9), X(12)]; dX = zeros(12,1);dX(1)=y;dX(2)=alf*cos(z)+a*x-c*x^3-b*y;dX(3)=omega;J=[0        1        0 ;           a-3*c*x^2   -b   -alf*sin(z);   0           0        0];dX(4:12) = J*Y;end


clear all;
clc;
% Z1=[ ];
% Z2=[ ];
% Z3=[ ];
global alf;
options=odeset('RelTol',1e-8,'AbsTol',1e-9);
for mm=1:100;
alf=(mm-1)*0.05
a0=0.1;b0=0.1
yinit = [a0,b0,1];
orthmatrix = eye(3);
IC(1:3) = yinit;
IC(4:12) = orthmatrix;
InitialTime = 0;                                           % 时间初始值
TimeStep = pi/180;                                         % 时间步长
wholetimes = 300*pi;                                       % 总的循环时间
UpdateStepNum =360;                                        % 每次演化的步数
iteratetimes = wholetimes/(UpdateStepNum*TimeStep);        % 演化的次数
DiscardItr=10;                                             % 抛弃的不稳定迭代次数
Iteration=UpdateStepNum*TimeStep;                          % 迭代次数
DiscardTime=DiscardItr*Iteration+InitialTime;              % 抛弃的时间段
mod = zeros(1,3);
d=3;                                                       % 维数
Sum=zeros(1,d);
n=0;                                                       % 总的迭代计数
k=0;                                                       % 对结果有作用的迭代计数
xData=[];
yData=[];
T1=InitialTime;
T2=T1+Iteration;
TSpan=[T1:TimeStep:T2];
        for i=1:iteratetimes
            n=n+1
            [t,X] =ode45('DF', TSpan, IC,options);
            [rX,cX]=size(X);   %提取行数和列数
            L=cX-d*d;      % 取积分得到的最后一个时刻的值
            for j=1:d
                m1=L+1+(j-1)*d;
                m2=m1+d-1;
                A(:,j)=(X(rX,m1:m2))';
            end
            y0 = ThreeGS(A);                      %正交归一化
            mod(1) = sqrt(y0(:,1)'*y0(:,1));   % 取三个向量的模
            mod(2) = sqrt(y0(:,2)'*y0(:,2));
            mod(3) = sqrt(y0(:,3)'*y0(:,3));
            y0(:,1) = y0(:,1)/mod(1);
            y0(:,2) = y0(:,2)/mod(2);
            y0(:,3) = y0(:,3)/mod(3);
             y(4:12) = y0';
            if (T2>DiscardTime)
                k=k+1;
                T=k*Iteration;
                TT=n*Iteration+InitialTime;

                %计算lyapunov指数
                Sum=Sum+log(abs(mod));
                lambda=Sum/T;
                %按降序排列指数
                Lambda=fliplr(sort(lambda));
                xData=[xData;TT];
                yData=[yData;lambda];
            end
            % 重新定义起始时刻
            ic=X(rX,1:L);
            T1=T2;
            T2=T2+Iteration;
            TSpan=[T1:TimeStep:T2];
            IC=[ic(:);y0(:)];
        end
        YY(mm,:)=Lambda'
end
本程序是从论坛下载，修改的，DF函数我还是不太明白另外结果好像不对，请赐教

Vickyvictoria

DF函数应该就是你所要计算的duffing方程

1713573225

请问你解决了吗我也在做这些可以想你请教下吗

christy

1713573225 发表于 2015-11-12 21:13
请问你解决了吗我也在做这些可以想你请教下吗

有问题建议直接在论坛提出来，大家一起讨论
因为楼主未必能够看到你发的内容