5.27

求解无约束优化问题： ；
（2）终止准则取 ；
（3）完成Newton法（牛顿法）的MATLAB编程、调试；
（4）选取几个与实验二中相同的初始点，并给出相关实验结果的对比及分析（从最优解、最优值、收敛速度（迭代次数）等方面进行比较）；
（5）按照模板撰写实验报告，要求规范整洁。
三、算法步骤、代码、及结果
1. 算法步骤
输入函数和参数，设置精度 epsilon， 转置初始点 x0 和变量 var_x 以确保它们是列向量。计算函数 f 的梯度 gradf 和 Hessian 矩阵 grad2f。初始化变量 grad_fxk（用于存储当前点的梯度值），设置初始迭代次数 k = 0，将当前点 xk 设为初始点 x0。
计算梯度和 Hessian 矩阵，迭代中如果梯度的范数 norm(double(grad_fxk)) 小于 epsilon，则停止迭代。最后返回优化后的解 x_optimization、优化后的函数值 f_optimization、总的迭代次数 k。多个初始点实验，对多个初始点 initial_points 进行实验，对每个初始点，调用 Multi_Newton_Method 并记录优化结果，保存结果到文件 experiment_results.mat。
2. 代码
function [x_optimization, f_optimization, k] = Multi_Newton_Method(f, x0, var_x, epsilon)

format long;
% 默认精度
if nargin == 3
epsilon = 1.0e-5;
end

% 转换初始点和变量
x0 = transpose(x0);
var_x = transpose(var_x);

% 计算梯度和 Hessian 矩阵
gradf = jacobian(f, var_x);
grad2f = jacobian(gradf, var_x);
grad_fxk = 1;
k = 0;
xk = x0;
while norm(double(grad_fxk)) > epsilon

% 计算梯度和 Hessian 在当前点的值
grad_fxk = subs(gradf, var_x, xk);
grad2_fxk = subs(grad2f, var_x, xk);

% 计算更新方向
pk = -inv(double(grad2_fxk)) * transpose(double(grad_fxk));
pk = double(pk);

% 更新变量
xk_next = xk + pk;
xk = xk_next;
k = k + 1;
end

% 返回优化结果
x_optimization = xk_next;
f_optimization = subs(f, var_x, x_optimization);
format short;
end

syms x1 x2 x3 x4
f = (x1 + 10*x2)^2 + 5*(x3 - x4)^2 + (x2 - 2*x3)^4 + 10*(x1 - x4)^4;
x0 = [1, 1, 1, 1]; % 选择实验二中的初始点
var_x = [x1, x2, x3, x4];
epsilon = 1e-6;

[x_optimization, f_optimization, k] = Multi_Newton_Method(f, x0, var_x, epsilon);

disp('优化后的解：');
disp(x_optimization);
disp('优化后的函数值：');
disp(f_optimization);
disp('迭代次数：');
disp(k);

% 选择多个初始点进行实验
initial_points = [
1, 1, 1, 1;
0, 0, 0, 0;
-1, -1, -1, -1;
2, 2, 2, 2
];

results = [];
for i = 1:size(initial_points, 1)
x0 = initial_points(i, :);
[x_optimization, f_optimization, k] = Multi_Newton_Method(f, x0, var_x, epsilon);
results = [results; x0, double(x_optimization'), double(f_optimization), k];
end

% 保存结果到文件
save('experiment_results.mat', 'results');

% 显示结果
disp('初始点\t\t优化后的解\t\t优化后的函数值\t\t迭代次数');
disp(results);