函数定义：
% 目标函数
function f = fun(x)
fm=0;
for i=1:499
fmi = (1-x(i))^2 + 100*(x(i+1)-x(i)^2)^2;
fm=fm+fmi;
end
f =fm;
end
% 梯度
function g = grad(x)
g = zeros(500,1);
g(1)=2*(x(1)-1)+400*x(1)*(x(1)^2-x(2));
for i=2:499
g(i)=2*(x(i)-1)+400*x(i)*(x(i)^2-x(i+1))+200*(x(i)-x(i-1)^2); end
g(500) = 200*(x(500)-x(499)^2);
end
% 二阶梯度
function g = grad2(x)
g = zeros(500,500);
g(1,1)=2+400*(3*x(1)^2-x(2));
g(1,2)=-400*x(1);
for i=3:500
g(1,i)=0;
end
for i=1:498
g(500,i)=0;
end
g(500,499)=-400*x(499);
g(500,500)=200;
for i=2:499
for j=1:500
if j==i-1
g(i,j)= -400*x(i-1);
elseif j==i
g(i,j)= 2+400*(3*x(i)^2-x(i+1))+200;
elseif j==i+1
g(i,j)= -400*x(i);
else g(i,j)=0;
end
end
end
end
1.最速下降法
function x_star = steepest(x0,eps)
gk = grad(x0);
res = norm(gk);
k = 0;
while res > eps && k<=50000
dk = -gk;
ak =1; f0 = fun(x0);
f1 = fun(x0+ak*dk);
slope = dot(gk,dk);
while f1 > f0 + 0.1*ak*slope
ak = ak/2;
xk = x0 + ak*dk;
f1 = fun(xk);
end
k = k+1;
x0 = xk;
gk = grad(xk);
res = norm(gk);
fprintf('--After the %d-th iter, the residual is %f\n',k,res);
end
x_star = xk;
end
运行结果：
>> steepest(zeros(500,1),0.0001)
--After the 1-th iter, the residual is 447.213595
--After the 2-th iter, the residual is 164.162920
……………………………
--After the 13587-th iter, the residual is 0.000111
--After the 13588-th iter, the residual is 0.000098
ans =
1.0000 （以下略，x各解均为1.0000）
2.牛顿法（不使用一维线搜索法，由于初始点离结果较远导致不收
敛，故采用初始值x均为0.8）
function x_star = newton(x0,eps)
gk = grad(x0);
bk = [grad2(x0)]^(-1);
res = norm(gk);
k = 0;
while res > eps && k<=1000
dk=-bk*gk;
xk = x0 + dk;
k = k+1;
x0 = xk;
gk = grad(xk);
bk = [grad2(xk)]^(-1);
res = norm(gk);
fprintf('--After the %d-th iter, the residual is %f\n',k,res); end
x_star = xk;
end
运行结果：
>> newton(0.8*ones(500,1),0.0001)
--After the 1-th iter, the residual is 95613.530371
--After the 2-th iter, the residual is 27958.861084
--After the 3-th iter, the residual is 8027.532086
--After the 4-th iter, the residual is 2197.358772
--After the 5-th iter, the residual is 524.525307
--After the 6-th iter, the residual is 82.202566
--After the 7-th iter, the residual is 4.168551
--After the 8-th iter, the residual is 0.010203
--After the 9-th iter, the residual is 0.002320
--After the 10-th iter, the residual is 0.000000
ans =
1.0000（以下略，x各解均为1.0000）
3.BFGS法
function x_star = bfgs(x0,eps)
g0 = grad(x0);
gk=g0;
res = norm(gk);
