本来想找一个现成的程序，结果找到的要么不能用，要么没有说明，索性自己编了一个。包括一个主程序和一个子程序。供大家分享

% Main program of steepest.

% 最速下降法主程序

%----------------------These lines only for ploting figures

%(前面几行是为了画图用的，可以不要)

clear

xx=-11:0.2:11; % 图中坐标的起始点

yy=-11:0.2:11;

[X,Y]=meshgrid(xx,yy);

for i=1:111

for j=1:111

y(1)=X(i,j);

y(2)=Y(i,j);

Z(i,j)=f_df(y); % z必须是网格化的数据

end

end

%mesh(X,Y,Z)

contour(X,Y,Z,50)

hold on

%-----------------------------------end for ploting figures

syms x1 x2 f lmda lm

x=[0;10]; % x的初值

xs=[x1; x2]; % 由于在求lmda过程中要用到符号函数的求导，故此定义符号运算的xs

[f0,g]=f_df(x); % 将x的初值代入，得到函数及其梯度的数值

[h, gs]=f_df(xs) % 将xs代入，得到符号表达式h

flag=999. % 用于控制循环

ii=0 % 用于循环计数

x_all=ones(1000,2); % 用于存储每次计算出来的x,如果循环步数超过1000，要增大数组

f_all=ones(1000,1); % 用于存储每次计算出来的f

g_all=ones(1000,2); % 用于存储每次计算出来的gradient

while flag>0.0001 % 当flag<0.01时退出循环

ii=ii+1

p=-g'; % 负梯度方向作为最速下降方向

xs(1)=x(1)+p(1)*lm;

xs(2)=x(2)+p(2)*lm; % 这里的lm实际上是lmda,但现在没有具体数值，只是一个符号

h1=compose(h,xs(1));

h2=compose(h1,xs(2)); % 这两步将f函数写成lmda的符号表达式

dh2=diff(h2,lm); % 对lmda符号函数求导

lmda=solve(dh2); % 用求导求极值的方法确定lmda

x_ini=double(x);

x_all(ii,:)=x_ini(:)';

f_all(ii)=double(f0);

g_all(ii,:)=double(g);

x=x+lmda*p; % 更新x

%--------------------------------------the nex two lines only for plot (这两行也是为画图用的)

x_end=double(x);

plot([x_ini(1),x_end(1)],[x_ini(2),x_end(2)],'-')

%-------------------

[f0,g]=f_df(x); % 更新f0和g

flag=sqrt(g(1)^2+g(2)^2); %已梯度的模作为控制循环的标记

flag=double(flag) %由于上一步是符号参数，不能直接放在while语句中控制，故此进行数值化

end

hold off

% subroution of steepest. Define a binary function and calculate it's gradient

% 最速下降法子程序，定义一个二元函数，并计算其梯度

function [f,df]=f_df(x);

f=5*x(1)^2+x(2)^2+4*x(1)*x(2)-14*x(1)-6*x(2); % for ellipse

df(1)=10*x(1)+4*x(2)-14;

df(2)=2*x(2)+4*x(1)-6;