|
- syms x1 x2;
- X = [x1;x2];
- fx = 2*x1^2 + 2*x1*x2 + 2*x2^2 - 4*x1 - 6*x2;
- x=[1;1];
- e=0.0001;
- minVal = newton(fx,X,x,e);
-
- function [ y ] = newton(fx,var,x,e,MAX)
- % 牛顿法求解函数极小点
- % 参数描述------------------------------
- % fx 符号表达式
- % var 符号变量列表 如:syms x1 x2;var= [x1;x2];
- % x 起始位置
- % e 精度控制
- % MAX 最大迭代次数控制
- % ------------------------------
-
- % Armijo步长
- step=2;% 初始步长设定
- r=0.5; % 取值范围(0,1),越大越快
- b=0.5; % 取值范围(0,1),越大越慢
-
- if nargin < 5
- MAX = 100;%设置默认最大迭代次数
- end
- precision = 5;%显示精度控制
- n=length(x);
-
- %开始循环迭代
- %direction存贮搜索方向
- %step 存贮步长
-
- H=eye(n);
- bfound = 0;
- for k=1:1:MAX
- [gx,direction] = getNextDirecrion(fx,var,x,H);
- disp('------------------------------');
- % fprintf('d[%d]=:',k);
- % disp( vpa(direction',precision) );
- if abs(subs(gx,var,x)) <= e
- y = x;
- bfound = 1;%得到结果时置为1
- break;
- else
- x0 = x+step*direction;
- while subs(fx,var,x0) > subs(fx,var,x) + b*step*r*subs(gx',var,x)*direction
- step = step*r;
- x0 = x+step*direction;
- end
- %打印求解过程
- fprintf('X[%d]=:',k);
- disp( vpa(x',precision) );
- fprintf('step(%d)=: ', k);
- disp( vpa(step,precision) );
- disp('------------------------------');
- x1 = x;
- x = x+step*direction;%计算下一个位置
- end
- H=getH(fx,var,x1,x,H);
- end
- if bfound == 1
- fprintf('X[%d]=:',k);
- disp( vpa(x',precision) );
- fprintf('step(%d)=: ', k);
- disp( vpa(step,precision) );
- disp('min value of:');
- fprintf("%.6f\n", vpa( subs(fx,var,y),precision) );
- end
- end
-
- %根据位置xk,获取搜索方向
- function [gx,direction] = getNextDirecrion(fx,var,xk,H)
- gx = gradient(fx,var); %计算梯度函数
- direction = -H^(-1)*subs(gx,var,xk);%计算搜索方向
- end
-
- %
- function [H] =getH(fx,var,x0,x1,H)
- gx = gradient(fx,var); %计算梯度函数
- s = x1 - x0;
- y = subs(gx,var,x1) - subs(gx,var,x0);
- H = H + y*y'/(y'*s) - H*s*s'*H/(s'*H*s);%BFGS
- % H = H + (y-H*s)*(y-H*s)'/((y-H*s)'*s);%秩1
- % H = H - H*y*y'*H/(y'*H*y) + s*s'/(s'*y); %DFP
- end
|