|
- 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 最大迭代次数控制
- % ------------------------------
-
- if nargin < 5
- MAX = 10;%设置默认最大迭代次数
- end
- precision = 3;%显示精度控制
- 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
- step = getNextStep(fx,var,x,direction);%计算步长
- if isempty(step)
- error('can not find a proper step.');
- end
- %打印求解过程
- fprintf('X[%d]=:',k);
- disp( vpa(x',precision) );
- fprintf('step(%d)=: ', k);
- disp( vpa(step,precision) );
- disp('------------------------------');
- x0 = x;
- x = x+step*direction;%计算下一个位置
- H=getH(fx,var,x0,x,H);
- end
- 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", 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
-
- %根据位置xk和方向dk,获取搜索步长step
- %注意符号表达式求导数的根时返回值转换为double类型
- function [step] =getNextStep(fx,var,xk,dk)
- syms lambda;
- phix = subs(fx,var,xk+lambda*dk);
- phix_diff = diff(phix);
- step = double(solve(phix_diff,'Real',true));%求取导函数的实数根
- 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
- end
|