5.29

共轭梯度法程序设计

function [x, k] = fr_conjugate_gradient()

    % 初始猜测

    x = [1; 2; 2; 2];

    

    % 计算初始梯度

    g = grad_f(x);

    

    % 初始搜索方向

    d = -g;

    

    % 终止准则的容差

    tol = 1e-6;

    

    % 最大迭代次数

    max_iter = 1000;

    

    for k = 1:max_iter

        % 线搜索找到最优步长 alpha

        alpha = line_search(@f, x, d);

        

        % 更新 x

        x = x + alpha * d;

        

        % 计算新的梯度

        g_new = grad_f(x);

        

        % 检查收敛性

        if norm(g_new) <= tol

            break;

        end

        

        % 计算 beta

        beta = (g_new' * g_new) / (g' * g);

        

        % 更新搜索方向

        d = -g_new + beta * d;

        

        % 更新梯度

        g = g_new;

    end

end

 

function y = f(x)

    % 计算目标函数值

    y = (x(1) + 10 * x(2))^2 + 5 * (x(3) - x(4))^2 + (x(2) - 2 * x(3))^4 + 10 * (x(1) - x(4))^4;

end

 

function g = grad_f(x)

    % 计算目标函数的梯度

    g = zeros(4,1);

    g(1) = 2 * (x(1) + 10 * x(2)) + 40 * (x(1) - x(4))^3;

    g(2) = 20 * (x(1) + 10 * x(2)) + 4 * (x(2) - 2 * x(3))^3;

    g(3) = 10 * (x(3) - x(4)) - 8 * (x(2) - 2 * x(3))^3;

    g(4) = -10 * (x(3) - x(4)) - 40 * (x(1) - x(4))^3;

end

 

function alpha = line_search(f, x, d)

    % 简单回溯线搜索

    alpha = 1;

    rho = 0.5;

    c = 1e-4;

    while f(x + alpha * d) > f(x) + c * alpha * (grad_f(x)' * d)

        alpha = alpha * rho;

    end

end