%非線性最小二乘法
clear;
fun=fun2;
syms x1 x2;
x0=[0;0];%設置初點
x=x0;
Dfx=[diff(fun,x1) diff(fun,x2)];
dfx=subs(Dfx,[x1 x2],x);
fx=subs(fun,[x1 x2],x);
e=1e-6;%設置精度精度
num=0;
d=-inv(dfx'*dfx)*dfx'*fx;
f=fun;
F=f(1)^2+f(2)^2;
alpha=1;
sigma=0.5;
rho=0.5;
DFx=[diff(F,x1);diff(F,x2)];
dFx=subs(DFx,[x1 x2],x);
while(subs(F,[x1 x2],x+alpha*d)>(subs(F,[x1 x2],x)+sigma*alpha*dFx'*d))
    alpha=rho*alpha;
end
num=num+1;
x_next=x+alpha*d;
while(sqrt((x_next-x)'*(x_next-x))>e)%開始用非線性最小二乘法迭代求最優點
    x=x_next;
    dfx=subs(Dfx,[x1 x2],x);
    fx=subs(fun,[x1 x2],x);
    d=-inv(dfx'*dfx)*dfx'*fx;
    alpha=1;
    sigma=0.5;
    rho=0.5;
    dFx=subs(DFx,[x1 x2],x);
    while(subs(F,[x1 x2],x+alpha*d)>(subs(F,[x1 x2],x)+sigma*alpha*dFx'*d))%開始用Armijo型線搜索方法求步長，并取σ=0.5，ρ=0.5
        alpha=rho*alpha;
    end
    x_next=x+alpha*d;
    num=num+1;
end
x=x_next;
Fx=subs(F,[x1 x2],x);
fprintf('\n非線性最小二乘法\n結果：\n優化精度為ε=%d\n初指點x0=[%d %d]\n最優點x=[ %d %d ] \nF(x)=%d \n||x(k+1)-x(k)||=%d \n迭代次數k=%d\n',e,x0(1),x0(2),x(1),x(2),Fx,sqrt((x_next-x)'*(x_next-x)),num);
[x,y]=meshgrid(-3:0.1:3,-3:0.1:3);
z=(x.^2-y+1).^2+(x-cos(3.1415926/2*y)).^2;
meshc(x,y,z);
xlabel('x1');
ylabel('x2');