工程数学实验三

function [x_opt, f_opt, iter] = newton_method()

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

    % 计算目标函数的梯度和Hessian矩阵
    grad_f = @(x) [400*x(1)*(x(1)^2 - x(2)) + 2*(x(1) - 1); -200*(x(1)^2 - x(2))];
    hessian_f = @(x) [1200*x(1)^2 - 400*x(2) + 2, -400*x(1); -400*x(1), 200];

    % 定义终止准则
    epsilon = 1e-5;

    % 设置初始点
    x0_list = [0, 0; -1, 1; 2, -2]; % 与上面的最速下降法使用相同的初始点

    for k = 1:size(x0_list, 1)
        % 初始化变量
        iter = 0;
        x_opt = x0_list(k, :)';
        f_opt = f(x_opt);
        grad_norm = norm(grad_f(x_opt));

        % 牛顿法迭代过程
        while grad_norm >= epsilon
            % 计算搜索方向
            p = -hessian_f(x_opt)\grad_f(x_opt);

            % 更新变量
            x_opt = x_opt + p;
            f_opt = f(x_opt);
            grad_norm = norm(grad_f(x_opt));

            % 更新迭代次数
            iter = iter + 1;
        end

        % 输出结果
        fprintf('初始点：%s\n', mat2str(x0_list(k, :)));
        fprintf('最优解：%s\n', mat2str(x_opt));
        fprintf('最优值：%f\n', f_opt);
        fprintf('迭代次数：%d\n', iter);
        disp('----------------------');
    end
end

% 调用牛顿法函数
[x_opt, f_opt, iter] = newton_method();