【最优化理论与算法】Newton法

syms x1 x2 x
import='输入原函数：';
f=input(import);
A=[0,0];                    
prompt='请输入初始点：';
A=input(prompt);             %初始点A
prompt='请输入允许误差：';
e=input(prompt);             %允许误差e
fx1=diff(f,'x1');            %计算梯度函数 [fx1,fx2]
fx2=diff(f,'x2');
fx1x1=diff(fx1,'x1');        % 计算黑塞矩阵函数 [fx1x1,fx1x2;fx2x1,fx2x2]
fx1x2=diff(fx1,'x2');
fx2x1=diff(fx2,'x1');
fx2x2=diff(fx2,'x2');
f1=subs(fx1,[x1,x2],[A(1),A(2)]);  %将A代入梯度函数 得到梯度df=[f1,f2]
f2=subs(fx2,[x1,x2],[A(1),A(2)]);
D=sqrt(f1^2+f2^2);   %计算误差D= || df ||
i=1;                 %迭代次数i
while D>e           %若误差D>精度e 继续迭代
    fprintf('第%d次迭代\n',i);
    i=i+1;
    f11=subs(fx1x1,[x1,x2],[A(1),A(2)]);  % 将A的值代入黑塞矩阵函数 得到矩阵的值
    f12=subs(fx1x2,[x1,x2],[A(1),A(2)]);
    f21=subs(fx2x1,[x1,x2],[A(1),A(2)]);
    f22=subs(fx2x2,[x1,x2],[A(1),A(2)]);
    df=[f1,f2];            %设梯度 df=[f1,f2];
    d2f=[f11,f12;f21,f22]; %黑塞矩阵 d2f=[f11,f12;f21,f22]
    d2f_n=inv(d2f);        %黑塞矩阵的逆d2f_n
    d=-df*d2f_n;           %计算搜索方向d
    fprintf('梯度为(%.6f,%.6f),搜索方向为(%.6f,%.6f),进行一维搜索\n',f1,f2,d(1),d(2));
    func=subs(f,[x1,x2],[A(1)+x*d(1),A(2)+x*d(2)]);  %得到沿d方向的一维搜索函数 func= min(A+x*d)
    k=One_dimensional_search(func);                  %用0.618法对函数一维搜索得到 k=λ
    A(1)=A(1)+k*d(1);                                %更新搜索初始点A
    A(2)=A(2)+k*d(2);
    fprintf('得到新的搜索点(%.6f,%.6f)\n',A(1),A(2));
    f1=subs(fx1,[x1,x2],[A(1),A(2)]);                %用新的A 更新梯度
    f2=subs(fx2,[x1,x2],[A(1),A(2)]);
    D=sqrt(f1^2+f2^2);                               %计算误差D= || df ||
    fprintf('误差为%.6f\n\n',D);
end
fprintf('误差小于允许误差，得到最优点(%.6f , %.6f)\n',A(1),A(2));
% (x1-1)^4+x2^2

function [y]=One_dimensional_search(func)  %% 使用0.618法的一维搜索函数，初始搜索区间由进退法获得
    [min_f,max_f]=fb(func,1);   % 进退法获得搜索区间 (min_f,max_f) 这里步长设置为1
    [u,v]=f_618(func,min_f,max_f,0.000001); % 0.618法获得符合精度的最优值区间 (u,v)  精度为10^-6
    y=(u+v)/2;   % 取平均值为函数最优点
    %%% 进退法，0.618法 函数代码块均在一维搜索函数内

    function [min_f,max_f]=fb(func,h0)   %%%%%%进退法函数代码块  func为函数，h0为步长
        h=h0;                     % 进退法step1  初始化
        x1=1;
        k=0;
        while 1
            x4 = x1 + h;          % 进退法step2  计算f(x4) f(x1)
            k = k + 1;
            f4 = subs(func,x4);
            f1 = subs(func,x1);
            if f4 < f1            % 进退法step4  f4<f1时执行 加大步长
                x2 = x1;
                x1 = x4;
                h = 2*h;
            else                  % 进退法step5  f4>f1时执行
                if k==1        % 进退法step6  改变步长方向
                    h = -h;       
                    x2 = x4;
                else           % 进退法step7  找到区间，退出函数
                    x3=x2;
                    x2=x1;
                    x1=x4;
                    min_f = min(x1, x3);
                    max_f = max(x1, x3);
                    break;  
                end
            end
        end
    end

    function [x1,x2]=f_618(func,a,b,L)  %%%%%     0.618法函数代码块 a,b为初始区间，由进退法获得，L为精度
        r=a+0.382*(b-a);           % step1 初始化
        s=a+0.618*(b-a);
        while 1    %重复以下过程，直至符合精度要求
            if b-a<L               %step2 达到精度要求，退出函数
                x1=a;
                x2=b;
                break;
            else                   %step2  未达到精度要求
                if(subs(func,r)>subs(func,s)) %step3 f(r)>f(s)时执行
                    a=r;
                    b=b;
                    r=s;
                    s=a+0.618*(b-a);
                else                          %step4 f(r)<f(s)时执行
                    a=a;
                    b=s;
                    s=r;
                    r=a+0.382*(b-a);
                end
            end
        end
    end
end