编写matlab程序使用最速下降算法、阻尼牛顿法、BFGS方法对Rosenbrock 函数、Beale 函数、Goldstein-Price 函数求解无约束优化问题。要求：正确求解出问题的最优解；画图展示三种方法各自的收敛路径；列表或画图比较算法的收敛速率；通过数值试验结果阐述三种方法各自的优缺点。

以下是matlab程序：

% Rosenbrock函数
f1 = @(x)(100 * (x(2) - x(1)^2)^2 + (1 - x(1))^2);
% Beale函数
f2 = @(x)((1.5 - x(1) + x(1)*x(2))^2 + (2.25 - x(1) + x(1)*x(2)^2)^2 + (2.625 - x(1) + x(1)*x(2)^3)^2);
% Goldstein-Price函数
f3 = @(x)((1 + (x(1) + x(2) + 1)^2 * (19 - 14*x(1) + 3*x(1)^2 - 14*x(2) + 6*x(1)*x(2) + 3*x(2)^2)) * (30 + (2*x(1) - 3*x(2))^2 * (18 - 32*x(1) + 12*x(1)^2 + 48*x(2) - 36*x(1)*x(2) + 27*x(2)^2)));

% 最速下降法
function [x, fval, iter] = gradient_descent(f, x0, alpha, eps)
    grad = @(x)([diff(f([x(1), x(2)] + [eps, 0]))/eps, diff(f([x(1), x(2)] + [0, eps]))/eps]);
    x = x0;
    iter = 0;
    while true
        iter = iter + 1;
        d = -grad(x);
        if norm(d) < eps
            break
        end
        x = x + alpha*d;
    end
    fval = f(x);
end

% 阻尼牛顿法
function [x, fval, iter] = newton_damp(f, x0, eps)
    grad = @(x)([diff(f([x(1), x(2)] + [eps, 0]))/eps, diff(f([x(1), x(2)] + [0, eps]))/eps]);
    hess = @(x)([diff(grad([x(1), x(2)] + [eps, 0]))/eps, diff(grad([x(1), x(2)] + [0, eps]))/eps]);
    x = x0;
    iter = 0;
    while true
        iter = iter + 1;
        d = -inv(hess(x))*grad(x)';
        if norm(d) < eps
            break
        end
        alpha = 1;
        while f(x+alpha*d') > f(x)
            alpha = alpha/2;
        end
        x = x + alpha*d';
    end
    fval = f(x);
end

% BFGS法
function [x, fval, iter] = BFGS(f, x0, eps)
    grad = @(x)([diff(f([x(1), x(2)] + [eps, 0]))/eps, diff(f([x(1), x(2)] + [0, eps]))/eps]);
    H = eye(2);
    x = x0;
    iter = 0;
    while true
        iter = iter + 1;
        d = -H*grad(x)';
        if norm(d) < eps
            break
        end
        alpha = 1;
        while f(x+alpha*d') > f(x)
            alpha = alpha/2;
        end
        s = alpha*d';
        x_old = x;
        x = x + s;
        y = grad(x) - grad(x_old);
        H = H + (y'*y)/(y'*s) - H*(s'*H*s)/(s'*H*s);
    end
    fval = f(x);
end

% 参数设置
alpha = 0.001;
eps = 1e-5;
x0 = [0, 0];

% Rosenbrock函数
[x_gd1, fval_gd1, iter_gd1] = gradient_descent(f1, x0, alpha, eps);
[x_nd1, fval_nd1, iter_nd1] = newton_damp(f1, x0, eps);
[x_bfgs1, fval_bfgs1, iter_bfgs1] = BFGS(f1, x0, eps);

% Beale函数
[x_gd2, fval_gd2, iter_gd2] = gradient_descent(f2, x0, alpha, eps);
[x_nd2, fval_nd2, iter_nd2] = newton_damp(f2, x0, eps);
[x_bfgs2, fval_bfgs2, iter_bfgs2] = BFGS(f2, x0, eps);

% Goldstein-Price函数
[x_gd3, fval_gd3, iter_gd3] = gradient_descent(f3, x0, alpha, eps);
[x_nd3, fval_nd3, iter_nd3] = newton_damp(f3, x0, eps);
[x_bfgs3, fval_bfgs3, iter_bfgs3] = BFGS(f3, x0, eps);

% 绘制收敛路径
[X, Y] = meshgrid(-2:0.05:2, -2:0.05:2);
Z1 = 100*(Y-X.^2).^2 + (1-X).^2;
Z2 = (1.5-X+X.*Y).^2 + (2.25-X+X.*Y.^2).^2 + (2.625-X+X.*Y.^3).^2;
Z3 = (1+(X+Y+1).^2.*(19-14*X+3*X.^2-14*Y+6*X.*Y+3*Y.^2)).*(30+(2*X-3*Y).^2.*(18-32*X+12*X.^2+48*Y-36*X.*Y+27*Y.^2));
figure(1)
contour(X, Y, Z1, 50)
hold on
plot(x_gd1(1), x_gd1(2), 'rx')
plot(x_nd1(1), x_nd1(2), 'bx')
plot(x_bfgs1(1), x_bfgs1(2), 'gx')
legend('Rosenbrock函数', '最速下降法', '阻尼牛顿法', 'BFGS法')
title('Rosenbrock函数收敛路径')

figure(2)
contour(X, Y, Z2, 50)
hold on
plot(x_gd2(1), x_gd2(2), 'rx')
plot(x_nd2(1), x_nd2(2), 'bx')
plot(x_bfgs2(1), x_bfgs2(2), 'gx')
legend('Beale函数', '最速下降法', '阻尼牛顿法', 'BFGS法')
title('Beale函数收敛路径')

figure(3)
contour(X, Y, Z3, 50)
hold on
plot(x_gd3(1), x_gd3(2), 'rx')
plot(x_nd3(1), x_nd3(2), 'bx')
plot(x_bfgs3(1), x_bfgs3(2), 'gx')
legend('Goldstein-Price函数', '最速下降法', '阻尼牛顿法', 'BFGS法')
title('Goldstein-Price函数收敛路径')

% 计算收敛速率
fval_gd1_list = zeros(1, iter_gd1);
fval_nd1_list = zeros(1, iter_nd1);
fval_bfgs1_list = zeros(1, iter_bfgs1);
for i = 1:iter_gd1
    [~, fval_gd1_list(i)] = gradient_descent(f1, x0, alpha, 10^(-i));
end
for i = 1:iter_nd1
    [~, fval_nd1_list(i)] = newton_damp(f1, x0, 10^(-i));
end
for i = 1:iter_bfgs1
    [~, fval_bfgs1_list(i)] = BFGS(f1, x0, 10^(-i));
end
figure(4)
semilogy(1:iter_gd1, abs(fval_gd1_list-fval_gd1), 'r')
hold on
semilogy(1:iter_nd1, abs(fval_nd1_list-fval_nd1), 'b')
semilogy(1:iter_bfgs1, abs(fval_bfgs1_list-fval_bfgs1), 'g')
legend('最速下降法', '阻尼牛顿法', 'BFGS法')
title('Rosenbrock函数收敛速率')

fval_gd2_list = zeros(1, iter_gd2);
fval_nd2_list = zeros(1, iter_nd2);
fval_bfgs2_list = zeros(1, iter_bfgs2);
for i = 1:iter_gd2
    [~, fval_gd2_list(i)] = gradient_descent(f2, x0, alpha, 10^(-i));
end
for i = 1:iter_nd2
    [~, fval_nd2_list(i)] = newton_damp(f2, x0, 10^(-i));
end
for i = 1:iter_bfgs2
    [~, fval_bfgs2_list(i)] = BFGS(f2, x0, 10^(-i));
end
figure(5)
semilogy(1:iter_gd2, abs(fval_gd2_list-fval_gd2), 'r')
hold on
semilogy(1:iter_nd2, abs(fval_nd2_list-fval_nd2), 'b')
semilogy(1:iter_bfgs2, abs(fval_bfgs2_list-fval_bfgs2), 'g')
legend('最速下降法', '阻尼牛顿法', 'BFGS法')
title('Beale函数收敛速率')

fval_gd3_list = zeros(1, iter_gd3);
fval_nd3_list = zeros(1, iter_nd3);
fval_bfgs3_list = zeros(1, iter_bfgs3);
for i = 1:iter_gd3
    [~, fval_gd3_list(i)] = gradient_descent(f3, x0, alpha, 10^(-i));
end
for i = 1:iter_nd3
    [~, fval_nd3_list(i)] = newton_damp(f3, x0, 10^(-i));
end
for i = 1:iter_bfgs3
    [~, fval_bfgs3_list(i)] = BFGS(f3, x0, 10^(-i));
end
figure(6)
semilogy(1:iter_gd3, abs(fval_gd3_list-fval_gd3), 'r')
hold on
semilogy(1:iter_nd3, abs(fval_nd3_list-fval_nd3), 'b')
semilogy(1:iter_bfgs3, abs(fval_bfgs3_list-fval_bfgs3), 'g')
legend('最速下降法', '阻尼牛顿法', 'BFGS法')
title('Goldstein-Price函数收敛速率')

三种方法各自的优缺点：

- 最速下降法：优点是算法简单易实现，缺点是收敛速度慢，在弯曲的目标函数上容易震荡。
- 阻尼牛顿法：优点是可以快速收敛到局部最优解，缺点是需要计算目标函数的海森矩阵，海森矩阵过大时计算代价高，且海森矩阵不一定正定，可以使用阻尼牛顿法解决。
- BFGS法：优点是收敛速度快，且不需要计算目标函数的海森矩阵，而是通过逼近海森矩阵来更新搜索方向，缺点是需要存储一个大小为$n^2$的矩阵，计算代价较高。