阻尼Newton法

发表于 更新于 分类于 阅读次数： 本文字数： 1.1k 阅读时长 ≈ 1 分钟

一般用来寻找二元函数最优解

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
function ZN_Newton(fun,x0,epsilon)
% 阻尼Newton法
%  fun = (1-x1)^2+2*(x2-x1^2)^2; x0 = [0;0]; epsilon = 0.1;
%  fun = (x1-1)^4+(x1-x2)^2; x0 = [0;0]; epsilon = 10^-6;
%  fun =  x1^2+x2^2-3*x1-x1*x2; x0 =[0;0];epsilon = 0.1;
%  fun =  (x1-2)^2+(x1-2*x2)^2; x0 =[0;3];epsilon = 0.4;
% clear;clc
% by -- 2021.12.9 -- cn.sgnxotsmicf.Lx
clc
disp('==================阻尼Newton法==================')
syms x1 x2
n = length(x0); % 解的长度
k = 0;    % 迭代次数
% 初始化值
grad_f = [diff(fun,x1);diff(fun,x2)];
Hesse = hessian(fun,[x1,x2]);
x = zeros(n,10); 
y = zeros(n,10);
index = 1;  % 索引
x(:,index) = x0;
% 进入迭代逻辑
while(true)
    syms lamuda real
    gk = subs(grad_f,[x1;x2],x(:,index));
    gk = double(gk);
    gk_Hesse = subs(Hesse,[x1;x2],x(:,index));
    y(:,index) = gk;
    if(norm(gk) > epsilon) % 精度判断
        dk = - gk_Hesse^-1 * gk;
        f1 = subs(fun,[x1;x2],x(:,index)+lamuda * dk);
        f1 = diff(f1);
        lamuda = solve(f1);      
                lamuda = max(lamuda); % 进入实数循环
                 if(lamuda ~= 0)
                      x(:,index+1) = x(:,index) + lamuda*dk;
                 else
                    break;
                 end      
    else
 % 迭代结束进行输出
     val = x(:,index);
     fprintf('阻尼Newton所求的函数为:%s\n迭代次数为;%d次\n',fun,k)
     str = ['最优解为:[',num2str(val'),']'];
     disp(str)
    break;
    end
    index = index + 1;
    k = k + 1;
end
end