henon.m

给出Henon的模型x(i+1)=1+y(i)-a*x(i)^2 y(i+1)=b*x(i)

%Henon模型的控制
clear all; tic;

t=[0:1:50];

%给出Henon的模型x(i+1)=1+y(i)-a*x(i)^2;y(i+1)=b*x(i);
x=zeros(1,length(t));y=zeros(1,length(t));
x0=0.5;y0=0.2;a=1.4;b=0.3;

%给出模拟退火算法的参数
T=10;La=10;K=0.75;beta=1.08;cita0=0.5;

%Henon模型的状态
S=2;
    for i=1:length(t)
        if i==1
            x(i)=x0;
            y(i)=y0;
        else
            if S==1
                x(i)=1+y0-a*x0^2;
                y(i)=b*x0;
            elseif S==2
                for j=1:La
                    cita1=rand(1);
                    X0=x0+cita0*(1+y0-a*x0^2-x0);
                    X1=x0+cita1*(1+y0-a*x0^2-x0);
                    dX0=abs((X0-x0)/x0);
                    dX1=abs((X1-x0)/x0);
                    dX=dX0-dX1;
                    if (dX>=0)|(exp(-abs(dX)/T)>rand(1))
                        cita0=cita1;
                    end
                end
                    x(i)=x0+cita0*(1+y0-a*x0^2-x0);
                    y(i)=b*x0;
                    T=K*T;La=beta*La;
            end
        end
        x0=x(i);
        y0=y(i);
    end

figure(1)
if S==1
    subplot(221);plot(t,x);title('Henon模型的混沌状态');subplot(222);plot(t,y);
elseif S==2
    subplot(223);plot(t,x);title('Henon模型的模拟退火控制');subplot(224);plot(t,y);
end
hold off
toc;