|
- syms x1 x2;
- X = [x1;x2];
- fx = 1.5*x1^2 - x1*x2 + 0.5*x2^2 - x1;
- x=[-2;4];
- e=0.01;
- 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;%显示精度控制
-
- %开始循环迭代
- %direction存贮搜索方向
- %step 存贮步长
-
- bfound = 0;
- for k=1:1:MAX
- [gx,direction] = getNextDirecrion(fx,var,x);
- 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('------------------------------');
- x = x+step*direction;%计算下一个位置
- 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\n", vpa( subs(fx,var,y),precision) );
- end
- end
-
- %根据位置xk,获取搜索方向
- function [gx,direction] = getNextDirecrion(fx,var,xk)
- gx = gradient(fx,var); %计算梯度函数
- G1=gradient(gx(1),var);
- G2=gradient(gx(2),var);
- H=[G1(1),G1(2);G2(1),G2(2)];%海森矩阵
- 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
|