5.29 matlab

% 定义目标函数
f = @(x) (x(1) + 10*2)^2 + 5*(x(3) - x(4))^2 + (x(2) - 2*x(3))^4 + 10*(x(1) - x(4))^4;

% 初始值和终止准则
x0_list = [-2, 2, -3, 3; -3, -1.5, 0.5, -1.5]; % 确保每个初始点有四个元素
tol = 1e-5;

% 梯度和海森矩阵函数（这里仅为示例，需要您根据实际情况计算正确的表达式）
grad_f = @(x) [2*(x(1) + 10*2); % 示例，实际梯度项需要根据f计算
                 4*(x(2) - 2*x(3))^3;
                 -10*(x(3) - x(4)) + 8*(x(2) - 2*x(3))^3;
                 10*4*(x(1) - x(4))^3 - 10*(x(3) - x(4))];

hes_f = @(x) [0, 0, 0, 0; 
               0, 12*(x(2) - 2*x(3))^2, -24*(x(2) - 2*x(3))^2, 0;
               0, -24*(x(2) - 2*x(3))^2, 40*(x(3) - x(4))^2 + 24*(x(2) - 2*x(3))^2, -40*(x(1) - x(4))^3 - 10;
               0, 0, -40*(x(1) - x(4))^3 - 10, 120*(x(1) - x(4))^2];

% 遍历每个初始点
for i = 1:size(x0_list, 1)
    x0 = x0_list(i, :);
    x = x0;
    iter = 0;
    grad_norm = inf;
    
    while grad_norm > tol
        iter = iter + 1;
        grad = grad_f(x);
        hes = hes_f(x);
        
        % 检查海森矩阵是否可逆
        if rank(hes) < length(x)
            warning('海森矩阵不可逆，可能遇到数值问题，退出迭代。');
            break;
        end
        
        % 求解牛顿方向
        d = -hes \ grad;
        
        % 更新x
        x_new = x + d;
        
        % 计算新的梯度范数
        grad_new = grad_f(x_new);
        grad_norm = norm(grad_new);
        
        % 检查是否收敛
        if isnan(grad_norm) || isinf(grad_norm)
            warning('梯度范数出现NaN或Inf，可能遇到数值问题，退出迭代。');
            break;
        end
        
        % 更新x
        x = x_new;
    end
    


fprintf('Initial point (%g, %g, %g, %g)\n', x0(1), x0(2), x0(3), x0(4));
    fprintf('Number of iterations: %d\n', iter);
    fprintf('Optimal point: (%g, %g, %g, %g)\n', x(1), x(2), x(3), x(4));
    fprintf('Optimal function value: %g\n', f(x));
    fprintf('\n');
end