无约束多维极值 阻尼牛顿法MATLAB代码程序

这个脚本包含两个主要的部分：一个是阻尼牛顿法（Min_ZuniNewton），另一个是黄金分割法（huangjin）。

阻尼牛顿法：

输入：目标函数f，初始点x0，允许的误差eps，最大迭代次数n。
输出：最小值点X和对应的函数值result。
算法步骤：首先，计算目标函数的梯度和海塞矩阵。然后，根据海塞矩阵是否为常数，选择合适的方式计算步长。接着，进行迭代，每次迭代使用牛顿法更新一步，然后使用黄金分割法确定步长。如果梯度的范数小于允许的误差，或者达到最大迭代次数，就停止迭代，返回当前点以及函数值。如果海塞矩阵是奇异的（条件数小于1e-6），就停止迭代，返回一个空数组和一个空结果。
黄金分割法：

输入：目标函数f，初始搜索区间[a, b]，初始步长h0，最大迭代次数n。
输出：最优步长result和最优步长对应的点x。
算法步骤：使用进退法找到初始的搜索区间，然后进行迭代。每次迭代，根据斐波那契序列生成两个新的候选点a和b，然后根据目标函数在这两个点的值，判断是否接受这两个点，缩小搜索区间。直到达到最大迭代次数，返回最优步长和对应的点。

%% 无约束多维极值  阻尼牛顿法
clc
clear
close
f=@(x1,x2) x1^2+2*x2^2-2*x1*x2-4*x1;
[X,result]=Min_ZuniNewton(f,[1 1],1e-8,100)
%阻尼牛顿法
function [X,result]=Min_ZuniNewton(f,x0,eps,n)
TiDu=gradient(sym(f),symvar(sym(f)));% 计算出梯度表达式
Haisai=jacobian(TiDu);
Var_Tidu=symvar(TiDu);
Var_Haisai=symvar(Haisai);
Var_Num_Tidu=length(Var_Tidu);
Var_Num_Haisai=length(Var_Haisai);
 
TiDu=matlabFunction(TiDu);
flag = 0;
if Var_Num_Haisai == 0  %也就是说海塞矩阵是常数
    Haisai=double((Haisai));
    flag=1;
end   
%求当前点梯度与海赛矩阵的逆
alfa_cal='f(';
f_cal='f(';
TiDu_cal='TiDu(';
Haisai_cal='Haisai(';
for k=1:length(x0)
    f_cal=[f_cal,'x0(',num2str(k),'),'];
    alfa_cal=[alfa_cal,'x1(',num2str(k),'),'];
    for j=1: Var_Num_Tidu
        if char(Var_Tidu(j)) == ['x',num2str(k)]
            TiDu_cal=[TiDu_cal,'x0(',num2str(k),'),'];
        end
    end
    
    for j=1:Var_Num_Haisai
        if char(Var_Haisai(j)) == ['x',num2str(k)]
            Haisai_cal=[Haisai_cal,'x0(',num2str(k),'),'];
        end
    end
end
Haisai_cal(end)=')';
TiDu_cal(end)=')';
f_cal(end)=')';
alfa_cal(end)=')';
 
switch flag
    case 0
        Haisai=matlabFunction(Haisai);
        dk='-eval(Haisai_cal)^(-1)*eval(TiDu_cal)';
    case 1
        dk='-Haisai^(-1)*eval(TiDu_cal)';
        Haisai_cal='Haisai';
end
 
i=1;
 
while i < n 
    if rcond(eval(Haisai_cal))<1e-6 % 判断海赛矩阵是否病态，实际上计算的是rcond(A),越小越病态
        disp('海赛矩阵奇异！');
        break;
    end
    syms alfa
    x1=x0(:)+alfa*eval(dk);%牛顿法迭代公式
    f_alfa=matlabFunction(eval(alfa_cal)); %计算得到包含alfa的表达式
    [~,alfa]=huangjin(f_alfa,0,0.1,100);%利用黄金分割法求得最优化步长因子αk
    x0=x0(:)+alfa*eval(dk);%拟牛顿法
    if norm(eval(TiDu_cal)) < eps
        X=x0;
        result=eval(f_cal);
        return;
    end
    i=i+1;
end
disp('无法收敛！');
X=[];
result=[];
end
%黄金分割法 
function [result,x]=huangjin(f,x0,h0,n)
tol=1e-6;
[a,b]=Min_jintui(f,x0,h0);%进退法寻找搜索区间
x1=min(a,b);
x2=max(a,b);
i=1;
while i < n
    %取中间值
    a=x1+0.382*(x2-x1);
    b=x1+0.618*(x2-x1);
    fa=f(a);
    fb=f(b);
    % 判断fa  fb大小，缩小区间
    if fa < fb
        x2=b;
    else
        x1=a;
    end
    if  abs(x1-x2) < tol
        result=f((x1+x2)/2);
        x=(x1+x2)/2;
        break;
    end
    i=i+1;
end
end
%进退法寻找搜索区间
function [a,b]=Min_jintui(f,x1,h0)
%f为函数句柄
%x1为初始点
%h0为步长
x2=x1+h0;
i=1;
if f(x2) < f(x1)
    x3=x2+h0;
    h=h0;
    while f(x3) < f(x2)
        x3=x2+h;
        h=i*h;
        i=i+1;
    end
    a=x1;
    b=x3;
else
    Temp=x1;
    x1=x2;
    x2=Temp;
    h=-h0;
    x3=x2+h;
    while f(x3) < f(x2)
        x3=x2+h;
        h=i*h;
        i=i+1;
    end
    a=x3;
    b=x1;
end
end