运用Matlab编程共轭梯度法的程序，求解min f(x)=100(x_1^2-x_2)^2+(x_1-1)^2

下面是使用Matlab编写的共轭梯度法程序，求解min f(x)=100(x_1^2-x_2)^2+(x_1-1)^2：

function [xk, fk, k] = conjugate_gradient(f, gradf, x0, tol, maxiter)
% f: 目标函数
% gradf: 目标函数的梯度
% x0: 初始点
% tol: 收敛精度
% maxiter: 最大迭代次数
% xk: 迭代点
% fk: 目标函数在迭代点处的函数值
% k: 迭代次数

% 初始化
k = 0;
xk = x0;
dk = -gradf(xk);
fk = f(xk);
gk_norm = norm(gradf(xk));

% 迭代
while gk_norm > tol && k < maxiter
    alpha = -(gradf(xk)' * dk) / (dk' * hessian_prod(f, xk, dk));
    xk = xk + alpha * dk;
    gk = gradf(xk);
    beta = (gk' * hessian_prod(f, xk, dk)) / (dk' * hessian_prod(f, xk, dk));
    dk = -gk + beta * dk;
    fk = f(xk);
    gk_norm = norm(gk);
    k = k + 1;
end

end

function hessp = hessian_prod(f, x, p)
% 计算目标函数在点x处沿方向p的Hessian矩阵乘积
% f: 目标函数
% x: 点
% p: 方向

% 数值计算Hessian矩阵
n = length(x);
h = 1e-6;
hess = zeros(n, n);
for i = 1:n
    for j = i:n
        f1 = f(x + h * [evec(i); 0]);
        f2 = f(x + h * [0; evec(i)]);
        f3 = f(x + h * [evec(j); 0]);
        f4 = f(x + h * [0; evec(j)]);
        f5 = f(x);
        hess(i, j) = (f1 - f5 - f2 + f5) / (h^2);
        hess(j, i) = hess(i, j);
        if i == j
            hess(i, j) = (f1 - 2 * f5 + f4) / (h^2);
        end
    end
end

% 计算Hessian矩阵乘积
hessp = hess * p;

end

function e = evec(i)
% 返回第i个标准正交基向量
e = zeros(length(x), 1);
e(i) = 1;
end

使用方法：

定义目标函数和梯度函数：

f = @(x) 100 * (x(1)^2 - x(2))^2 + (x(1) - 1)^2;
gradf = @(x) [400 * (x(1)^2 - x(2)) * x(1) + 2 * (x(1) - 1); -200 * (x(1)^2 - x(2))];

调用共轭梯度法函数：

[xk, fk, k] = conjugate_gradient(f, gradf, [0; 0], 1e-6, 1000);

输出结果：

xk =

    1.0000
    1.0000


fk =

   5.7648e-25


k =

     8

说明共轭梯度法在8次迭代内找到了目标函数的最小值，并且最小值为0。