MATLAB实现牛顿法（Newton‘s method）寻优

这是一个用MATLAB实现的牛顿法（Newton's method）的函数。牛顿法是一种在实数域和复数域上近似求解方程的方法。它使用函数f的泰勒级数的前几项来寻找方程f(x)=0的根。

在这个函数中，输入参数包括：

f：要求解的函数。
x0：初始的猜测值。
eps：精度要求。
n：最大迭代次数。
函数的输出是求解的函数f在x0处的最小值及其对应的x值。

下面是代码的详细解释：

定义函数f，输入为两个变量x1和x2。
定义Min_Newton函数，输入为要求解的函数f，初始猜测值x0，精度eps和最大迭代次数n。
在Min_Newton函数中，首先计算出函数f的梯度表达式TiDu和海塞矩阵表达式Haisai。
然后将梯度表达式和海塞矩阵表达式转换为MATLAB函数。
根据海塞矩阵表达式中是否有参数来选择合适的语句。
在while循环中，进行最大迭代次数n次的迭代。在每次迭代中，根据牛顿法的原理，计算新的猜测值x0，如果梯度的模小于eps，就认为找到了最小值，退出循环。
如果超过了最大迭代次数还没有找到满足精度要求的结果，就认为无法收敛。

这个函数实现的是基本的牛顿法，对于一些特殊情况如海塞矩阵是奇异矩阵或者病态矩阵等没有进行处理。在实际使用中，可能需要根据具体情况进行修改和优化。

%牛顿法
f=@(x1,x2) x1^2+2*x2^2-2*x1*x2-4*x1;
[X,result]=Min_Newton(f,[1 1],0.1,100)
function [X,result]=Min_Newton(f,x0,eps,n)
 
TiDu=gradient(sym(f),symvar(sym(f)));% 计算出梯度表达式
Haisai=jacobian(TiDu,symvar(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
%求当前点梯度与海赛矩阵的逆
f_cal='f(';
TiDu_cal='TiDu(';
Haisai_cal='Haisai(';
for k=1:length(x0) %求得初始变量的x0的元素个数
    f_cal=[f_cal,'x0(',num2str(k),'),'];%组装f_cal=f(x0(k))求得该点函数值
  
    for j=1: Var_Num_Tidu %求得梯度中的元素个数
        if char(Var_Tidu(j)) == ['x',num2str(k)] 
            TiDu_cal=[TiDu_cal,'x0(',num2str(k),'),'];%组装TiDu_cal=TiDu_cal(x0(k)求得该点梯度值
        end
    end
    
    for j=1:Var_Num_Haisai
        if char(Var_Haisai(j)) == ['x',num2str(k)]
            Haisai_cal=[Haisai_cal,'x0(',num2str(k),'),'];%组装Haisai_cal=Haisai_cal(x0(k)求得该点海塞矩阵的值
        end
    end
end
Haisai_cal(end)=')';  %完成海塞矩阵封装
TiDu_cal(end)=')';%完成梯度封装 
f_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 %设置最大迭代次数n
    if rcond(eval(Haisai_cal)) < 1e-6 %计算海塞矩阵的条件数 条件数越大，逆阵结果越不稳定
        disp('海赛矩阵病态！'); %条件数超出范围，示为病态矩阵
        break;
    end
    x0=x0(:)+eval(dk);   %eval函数将字符串转换为指令
    if norm(eval(TiDu_cal)) < eps 
        X=x0;
        result=eval(f_cal); 
        return;
    end
    i=i+1;
end
disp('无法收敛！'); %超出迭代范围
X=[];
result=[];
end