6.5

实验三：Newton法程序设计

一、实验目的

掌握Hesse矩阵的计算方法和Newton法的基本思想及其迭代步骤；学会运用MATLAB编程实现常用优化算法；能够正确处理实验数据和分析实验结果及调试程序。

二、实验内容

（1）求解无约束优化问题：f(x)=(x1+10x2)^2+5(x3-x4)^2+(x2-2x3)^4+10(x1-x4)^4；

（2）终止准则取||眩f(x^k)|| 10^-6；

（3）完成Newton法（牛顿法）的MATLAB编程、调试；

（4）选取几个与实验二中相同的初始点，并给出相关实验结果的对比及分析（从最优解、最优值、收敛速度（迭代次数）等方面进行比较）；

（5）按照模板撰写实验报告，要求规范整洁。

三、算法步骤、代码、及结果

function newton_unconstrained_optimization

% 定义目标函数 f 及其梯度 grad_f 和 Hessian 矩阵 hessian_f

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

grad_f = @(x) [

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

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

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

-10*(x(3) - x(4)) - 40*(x(1) - x(4))^3

];

hessian_f = @(x) [

2 + 120*(x(1) - x(4))^2, 20, 0, -120*(x(1) - x(4))^2;

20, 200 + 12*(x(2) - 2*x(3))^2, -24*(x(2) - 2*x(3))^2, 0;

0, -24*(x(2) - 2*x(3))^2, 10 + 48*(x(2) - 2*x(3))^2, -10;

-120*(x(1) - x(4))^2, 0, -10, 10 + 120*(x(1) - x(4))^2

];

% 参数设置

epsilon = 1e-6; % 终止准则的阈值

max_iter = 1000; % 最大迭代次数

initial_points = [

1, 1, 1, 1;

-1, 2, -3, 4;

0.5, -0.5, 0.5, -0.5

]; % 多个初始点

% 记录每个初始点的结果

results = cell(size(initial_points, 1), 1);

for i = 1:size(initial_points, 1)

x = initial_points(i, :)';

iter = 0;

f_values = [];

while norm(grad_f(x)) > epsilon && iter < max_iter

iter = iter + 1;

hessian = hessian_f(x);

gradient = grad_f(x);

% 求解线性系统 H * p = -grad 以得到搜索方向 p

p = -hessian \ gradient;

alpha = line_search(f, grad_f, x, p); % 线搜索确定步长

x = x + alpha * p;

f_values(end+1) = f(x); % 记录函数值

end

results{i} = struct('x', x, 'f_value', f(x), 'iterations', iter, 'f_values', f_values);

end

% 打印结果

for i = 1:length(results)

fprintf('Initial point: (%.2f, %.2f, %.2f, %.2f)\n', initial_points(i, 1), initial_points(i, 2), initial_points(i, 3), initial_points(i, 4));

fprintf('Optimal solution: (%.6f, %.6f, %.6f, %.6f)\n', results{i}.x);

fprintf('Optimal function value: %.6f\n', results{i}.f_value);

fprintf('Iterations: %d\n\n', results{i}.iterations);

end

% 画出最优值随迭代次数变化的曲线图

for i = 1:length(results)

figure;

plot(results{i}.f_values);

title(sprintf('Initial point: (%.2f, %.2f, %.2f, %.2f)', initial_points(i, 1), initial_points(i, 2), initial_points(i, 3), initial_points(i, 4)));

xlabel('Iterations');

ylabel('Function Value');

end

end

% 线搜索函数

function alpha = line_search(f, grad_f, x, p)

alpha = 1;

c = 1e-4;

rho = 0.9;

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

alpha = rho * alpha;

end

end