1.题目

设初始点为[x1 x2]=[-2.5 4.5] ，ε≤0.01，求目标函数的极小值。

2.matlab代码

2.1主函数

syms x1 x2 s; %声明符号变量

f1 =  x1^4 - 2*x1^2*x2 - 2*x1*x2 + x1^2 + 2*x2^2 + 4.5*x1 - 4*x2 + 4;%设定目标函数

k=steepest_descent(f1,x1,x2,s,[-2.5,4.25],10^(-2));  %设定起始点[x1 x2]=[-2.5,4.25]和精度10^(-2)

result_string=sprintf('在 %d 次迭代后求出极小值\n',k);%在迭代多少次之后求出极小值
disp(result_string);

2.2调用函数

function   k = steepest_descent(f,x1,x2,s,start_point,thereshold) 
    
    k = 0;%迭代次数赋值初始化
    
    grad_f = [diff(f,x1) diff(f,x2)]; %计算f的梯度
    
    delta = subs(grad_f,[x1,x2],[start_point(1),start_point(2)]);
    %计算起点的梯度
    
    step=1; %设置初始步长为1
    
    current_point = start_point;%起点值赋给当前点
    
    %最速下降法的主循环，判断条件为：梯度的模与所给精度值进行比较
    while norm(delta) > thereshold  
        
        k = k + 1;%迭代次数+1
        
        %一维探索 求最优步长（此时方向已知，步长s为变量）
        x_next = [current_point(1),current_point(2)] - s* delta/norm(delta);% 计算x（k+1）点，其中步长s为变量 
        f_val = subs(f,[x1,x2],[x_next(1),x_next(2)]);% 将x值带入目标函数中
        step = abs(double(solve(diff(f_val,s)))); % 对s求一阶导，并加绝对值符号，得到最优步长的绝对值
        step = step(1);%更新步长
        
        %计算x（k+1）点
        current_point = double([current_point(1),current_point(2)] - step * delta/norm(delta));
        
        %计算x（k+1）点的梯度值
        delta = subs(grad_f,[x1,x2],[current_point(1),current_point(2)]);
        
        %计算函数值
        f_value = double(subs(f,[x1,x2],[current_point(1),current_point(2)]));
        
        %输出迭代计算过程
        result_string=sprintf('k=%d, x1=%.6f, x2=%.6f, step=%.6f f(x1,x2)=%.6f',...
        k,current_point(1),current_point(2),step,f_value);
        disp(result_string);
        
    end
end

（使用代码时，只需要把主函数中目标函数 ；起始点； 精度，重新设置成你所需要的就可运行。）

2.3运行结果

3.分析

将收敛条件由ε≤0.01改为ε≤10^(-6)，运行结果如下：

可以看出，利用最速下降法求函数极小值时，在最初几步迭代中函数值下降很快，但当函数值越接近理论极小值时，函数值下降的越慢，同时，越接近理论极小值时，步长也越小，因此最速下降法的收敛速度并不快，这是因为函数在x（k）点处的负梯度方向为其最速下降方向仅仅是针对该点处而言，，一旦离开该点，原先方向就不是最速下降方向了。

最速下降方向只是针对当前的计算点而言，并非是全局的最速下降方向

最速下降法优点：对初始点选择要求低，远离极值点时收敛速度快
     缺点：越逼近理论极小值，收敛速度越慢

matlab实现最速下降法
定义：沿负梯度方向进行搜索的算法（负梯度方向为最速下降方向）

算法步骤：

步0：选取初始点x0，容许误差是e=[0~1],令k=1

步1：计算目标函数的梯度.若 ||gk||<=e，即达到误差要求,立即停止计算，并输出xk作为近似最优解。

步2:取搜索方向为dk=-gk（即负梯度方向）。

步3：利用线搜索技术确定步长k（这里采用Armijo准则来求步长） 

步长为k=^mk       是给定的，所以要求出mk

Amrijo准则就是

（1）给定（0~1），（0，0.5），令m=0

（2）若不等式

f(xk+^m*dk)<=f(xk)+*^m*gk'*dk

成立，则令mk=m，Xk+1=xk+^m*dk.停止运算，输出mk得到步长

（3）若不满足上述不等式，则令m=m+1，然后回到第二步。

步4：确定步长后，令Xk+1=Xk+k*dk，k=k+1，转步骤1.

matlab 具体代码如下：

1.主函数

clear all
clc
%利用grad函数求解 minif(x)=100*(x1^2-x2)^2+(x1-1)^2
%此时还要建立两个函数，一个目标函数fun,一个梯度gfun
x0=[-1.2 1]';
[x,val,k]=grad('fun','gfun',x0);
disp(['最优解：x = '])
disp(x)
disp(['此时: f(x) = ',num2str(val)]) 

 2.最速下降法

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

 3.目标函数

function f= fun(x)
%目标函数
f=100*(x(1)^2-x(2))^2+(x(1)-1)^2;
end

4.目标函数的梯度

function  g=gfun(x)
%目标函数的梯度
g=[400*x(1)*(x(1)^2-x(2))+2*(x(1)-1),-200*(x(1)^2-x(2))]';
end

5.运行结果