%%初始化
clear all;clc;
%構造函數
[x,y] = meshgrid([-4:0.5:4]);
z=(sin(x.*y)+x.^2+y.^2);
%初始化各個參數
k=1.3806505e-23 ;alpha=0.95;
%初始溫度
T=1000;
%給一個隨機位置
xx=(rand-1/2)*8;yy=(rand-1/2)*8;
%記錄坐標
xxx=[xx];yyy=[yy];
%%退火算法
for i=1:10000
    deltax=(rand-0.5)*0.04;deltay=(rand-0.5)*0.04;
    deltaz=sin((xx+deltax)*(yy+deltay))+(xx+deltax)^2+(yy+deltay)^2 ...
          -(sin(xx*yy)+xx^2+yy^2);
    P=exp(-deltaz/(k*T));
    if sign(deltaz) == -1
        xx=xx+deltax;
        yy=yy+deltay;        
        else if sign(rand-P) == -1
                xx=xx+deltax;
                yy=yy+deltay;
             end
    end
    T=alpha*T;
    %%xxx=[xxx xx];
    %%yyy=[yyy yy];
    if i/100 == fix(i/100)
        xxx=[xxx xx];
        yyy=[yyy yy];
    end
end          
subplot(2,2,1);surf(x,y,z);
subplot(2,2,2);plot3(x,y,z);hold on;plot3(xxx,yyy,sin(xxx.*yyy)+xxx.^2+yyy.^2,'-*');
subplot(2,2,3);plot3(xxx,yyy,sin(xxx.*yyy)+xxx.^2+yyy.^2,'-*');
subplot(2,2,4);plot(xxx,yyy);