matlab中牛顿法编程

Posted 2023-05-15

要计算一个很复杂的矩阵H（含有未知数x），打算通过det（H）=0之后，求出x的数值。因为矩阵中含有很多特殊函数（比如贝塞尔函数），所以用solve求解的时候，多数会求解不出来。所以改用了牛顿法。
但是在最近的计算中发现，那些能用solve求解出的答案和用牛顿法求解出的答案差好几倍。所以不知道是牛顿法的问题，还是因为矩阵太复杂用solve计算已经不可靠的问题？
所以恳求大家帮忙分析一下是为什么，并编写一个能求出正确解的程序（最好是用牛顿法）。多谢大家！

这样吧，你不要直接通过solve求解啊，可以画一个图像啊，以x的值为横坐标，行列式值为纵坐标，或者纵坐标取一下对数（如果变化范围太大），然后可以得出大致得到解的范围；最后可以通过其他方法计算更加精确的解。

希望能帮到你，我算程序就是这样实现的 参考技术A 我也期待答案。

第三章 最速下降法和牛顿法

内容来自马昌凤编著的《最优化方法及其Matlab程序设计》，文章仅为个人的学习笔记，感兴趣的朋友详见原书。
本章讨论无约束优化问题$minf(x)$的最速下降法和牛顿法及其改进算法，前者简单而古老，虽不再具有实用性，却是研究其他无约束优化算法的基础；后者也是一种经典的无约束优化算法，具有收敛速度快和自适应性等优点。

最速下降方法及其Matlab实现

算法

程序

function [x,val,k]=grad(fun,gfun,x0)
% 功能: 用最速下降法求解无约束问题:  min f(x)
%输入:  x0是初始点, fun, gfun分别是目标函数和梯度
%输出:  x, val分别是近似最优点和最优值,  k是迭代次数.
maxk=5000;   %最大迭代次数
rho=0.5;sigma=0.4;
k=0;  epsilon=1e-5;
while(k<maxk)
    g=feval(gfun,x0);  %计算梯度
    d=-g;    %计算搜索方向（精确线搜索）
    if(norm(d)<epsilon), break; end
    m=0; mk=0;
    while(m<20)   %Armijo搜索
        if(feval(fun,x0+rho^m*d)<feval(fun,x0)+sigma*rho^m*g'*d)
            mk=m; break;
        end
        m=m+1;
    end
    x0=x0+rho^mk*d;
    k=k+1;
end
x=x0;
val=feval(fun,x0); 

注：这里的$fun$与$gfun$分指目标函数与梯度，需在grad运行前准备好。（Armijo搜索部分可详见第二章相关内容）

示例1

示例2

牛顿法及其Matlab实现

基本思想：用迭代点$x_k$处的一阶导数和二阶导数对目标函数进行二次函数近似（泰勒多项式），然后把二次模型的极小点作为新的迭代点，并不断重复这一过程，直至求得满足精度的近似极小点。

算法

牛顿法最突出的优点是收敛速度快，具有局部二阶收敛性。但是对初始点要求足够“靠近”极小点，否则，有可能导致算法不收敛。由于实际问题的精度极小点一般是未知的，为克服这一困难，故引入线搜索技术以得到大范围收敛的算法，即下面的阻尼牛顿法。（此处基于Armijo准则）

阻尼牛顿法

算法

程序

function [x,val,k]=dampnm(fun,gfun, Hess,x0)
%功能: 用阻尼牛顿法求解无约束问题:  min f(x)
%输入: x0是初始点, fun, gfun, Hess 分别是求
%         目标函数,梯度,Hesse 阵的函数
%输出:  x, val分别是近似最优点和最优值,  k是迭代次数.
maxk=100;   %给出最大迭代次数
rho=0.55;sigma=0.4;
k=0;  epsilon=1e-5;
while(k<maxk)
    gk=feval(gfun,x0); %计算梯度
    Gk=feval(Hess,x0);  %计算Hesse阵
    dk=-Gk\\gk; %解方程组Gk*dk=-gk, 计算搜索方向
    if(norm(gk)<epsilon), break; end  %检验终止准则
    m=0; mk=0;
    while(m<20)   % 用Armijo搜索求步长 
        if(feval(fun,x0+rho^m*dk)<feval(fun,x0)+sigma*rho^m*gk'*dk)
            mk=m; break;
        end
        m=m+1;
    end
    x0=x0+rho^mk*dk;
    k=k+1;
end
x=x0;
val=feval(fun,x); 
%gval=norm(gfun(x));

此处除需提前建立目标函数和梯度的M文件外，还需建立Hess矩阵的求解文件

function He=Hess(x)
n=length(x)
He=zeros(n,n);
He=[此处填函数的二阶求导结果]

修正牛顿法及其Matlab实现

牛顿法具有不低于二阶的收敛速度，但该算法要求目标函数的Hess矩阵在每个迭代点$x_k$处是正定的，否则，难以保证牛顿方向$d_k=-G_k^{-1}{g}_k$是$f$在$x_k$处的下降方向。为克服这一缺陷，故对其进行修正：将牛顿法与最速下降法结合起来，构造**“牛顿-最速下降算法”**。其基本思想为：当Hess正定时，采用牛顿方向作为搜索方向；否则，采用负梯度方向作为搜索方向。

算法

程序

function [x,val,k]=revisenm(fun,gfun,Hess,x0)
%功能: 用修正牛顿法求解无约束问题:  min f(x)
%输入: x0是初始点, fun, gfun, Hess 分别是求
%         目标函数,梯度,Hesse 阵的函数
%输出:  x, val分别是近似最优点和最优值,  k是迭代次数.
n=length(x0); maxk=150;
rho=0.55;sigma=0.4; tau=0.0;
k=0;  epsilon=1e-5;
while(k<maxk)
    gk=feval(gfun,x0); % 计算梯度
    muk=norm(gk)^(1+tau);
    Gk=feval(Hess,x0);  % 计算Hesse阵
    Ak=Gk+muk*eye(n);
    dk=-Ak\\gk; %解方程组Gk*dk=-gk, 计算搜索方向
    if(norm(gk)<epsilon), break; end  %检验终止准则
    m=0; mk=0;
    while(m<20)   %用Armijo搜索求步长 
        if(feval(fun,x0+rho^m*dk)<feval(fun,x0)+sigma*rho^m*gk'*dk)
            mk=m; break;
        end
        m=m+1;
    end
    x0=x0+rho^mk*dk;
    k=k+1;
end
x=x0;
val=feval(fun,x); 
%gval=norm(gfun(x));