MATLAB的牛顿法求多元函数的极值程序加实例

　　今天主要是讲解MATLAB的牛顿法求多元函数的极值程序加实例。
　　实例1
　　求f(x,y)= sin(x^2+y^2)*exp(-0.1*(x^2+y^2+x*y+2*x))，在-2<=x<=2，-2<=y<=2上的极值点和极值。
　　主程序 clc; clear all; close all; syms x y;%定义函数变量 x y f = sin(x^2+y^2)*exp(-0.1*(x^2+y^2+x*y+2*x)); x0 = [1 1];%初始点 x0(1,1) [x_best,f_best] = Newton(f,x0,[x y]); x_best f_best = vpa(f_best) x = -2:0.01:2; y = x; [X,Y] = meshgrid(x,y); F = sin(X.^2+Y.^2)*exp(-0.1*(X.^2+Y.^2+X.*Y+2.*X)); figure; mesh(X,Y,F); xlabel(＂x＂); ylabel(＂y＂); zlabel(＂z＂);
　　牛顿法函数 function [x_best,f_best] = Newton(f,x0,x,epsilon) %% 牛顿法求解函数的最小值(极小值) %% 输入 %   f：目标函数 %   x0：初始点 %   x：自变量向量 %   epsilon：精度 %% 输出 %   x_bes：目标函数取最小值时的自变量值 %   f_best：目标函数的最小值 format long;%改变数据显示格式 if nargin == 3  %默认的精度     epsilon = 1.0e-6; end x0 = transpose(x0);%transpose函数的功能是转置向量或矩阵 x = transpose(x);%transpose函数的功能是转置向量或矩阵 g1f = jacobian(f,x);% jacobian求解向量函数的雅可比矩阵式  g2f = jacobian(g1f,x);% jacobian求解向量函数的雅可比矩阵式  % 参数初始化 grad_fxk = 1; k = 0; xk = x0;   while norm(grad_fxk) > epsilon  %   计算矩阵 (向量) X的2-范数     grad_fxk  = subs(g1f,x,xk);%   计算矩阵 (向量) 雅可比矩阵式在xk处的值     grad2_fxk = subs(g2f,x,xk);     pk = -inv(grad2_fxk)*transpose(grad_fxk);  % 步长     pk = double(pk);%转化为双精度浮点类型     xk_next = xk + pk; %       xk = xk_next;     k = k + 1;     f_1 = subs(f,x,xk);%计算函数值     %输出迭代结果     fprintf(＂迭代次数：%d  误差：%.20f 极值点：(x,y) = (%f,%f) 极值：f(x,y) = %.20f ＂,k,vpa(norm(grad_fxk)),xk(1),xk(2),vpa(f_1)); end %输出极值点和极值 x_best = xk_next; f_best = subs(f,x,x_best); end
　　运行结果 迭代次数：1  误差：1.02885710610701086587 极值点：(x,y) = (0.669084,0.966374) 极值：f(x,y) = 0.70142228466448164337 迭代次数：2  误差：0.14448082736806977522 极值点：(x,y) = (1.195944,0.595077) 极值：f(x,y) = 0.59942448686119498280 迭代次数：3  误差：0.67873695620313101440 极值点：(x,y) = (1.032695,0.554239) 极值：f(x,y) = 0.65658602325338621952 迭代次数：4  误差：0.03278835230868389766 极值点：(x,y) = (1.077563,0.457762) 极值：f(x,y) = 0.65569150404015985600 迭代次数：5  误差：0.01819636638003245543 极值点：(x,y) = (1.069052,0.464828) 极值：f(x,y) = 0.65572832791085189363 迭代次数：6  误差：0.00027874333536557117 极值点：(x,y) = (1.069330,0.464057) 极值：f(x,y) = 0.65572826847418552720 迭代次数：7  误差：0.00000108627104183494 极值点：(x,y) = (1.069329,0.464058) 极值：f(x,y) = 0.65572826847430654151 迭代次数：8  误差：0.00000000000108544724 极值点：(x,y) = (1.069329,0.464058) 极值：f(x,y) = 0.65572826847430654151   x_best =      1.069329230413560    0.464057718471801     f_best =   0.65572826847430659287489727298377
　　实例2
　　求f(x,y)= 4*(x-y)-x^2-y^2，在-2<=x<=2，-2<=y<=2上的极值点和极值。
　　主程序 clc; clear all; close all; syms x y;%定义函数变量 x y fx = 4*(x-y)-x^2-y^2;%定义二元变量函数 x0 = [1 1];%初始点 x0(1,1) [x_best,f_best] = Newton(fx,x0,[x y]); x_best f_best = vpa(f_best) x = -2:0.1:2; y = x; [X,Y] = meshgrid(x,y); F =  4.*(X-Y)-X.^2-Y.^2; figure; mesh(X,Y,F); xlabel(＂x＂); ylabel(＂y＂); zlabel(＂z＂);
　　运行结果 迭代次数：1  误差：6.32455532033675904557 极值点：(x,y) = (2.000000,-2.000000) 极值：f(x,y) = 8.00000000000000000000 迭代次数：2  误差：0.00000000000000000000 极值点：(x,y) = (2.000000,-2.000000) 极值：f(x,y) = 8.00000000000000000000   x_best =        2     -2     f_best =   8.0
　　实例3
　　求f(x,y)= (1-x)^2+100*(y-x^2)^2，在-2<=x<=2，-2<=y<=2上的极值点和极值。
　　主程序 clc; clear all; close all; syms x y;%定义函数变量 x y f = (1-x)^2+100*(y-x^2)^2; x0 = [0 0];%初始点 x0(1,1) [x_best,f_best] = Newton(f,x0,[x y]); x_best f_best = vpa(f_best) x = -2:0.1:2; y = x; [X,Y] = meshgrid(x,y); F = (1-X).^2+100.*(Y-X.^2).^2; figure; mesh(X,Y,F); xlabel(＂x＂); ylabel(＂y＂); zlabel(＂z＂);
　　运行结果 迭代次数：1  误差：2.00000000000000000000 极值点：(x,y) = (1.000000,0.000000) 极值：f(x,y) = 100.00000000000000000000 迭代次数：2  误差：447.21359549995793258859 极值点：(x,y) = (1.000000,1.000000) 极值：f(x,y) = 0.00000000000000000000 迭代次数：3  误差：0.00000000000000000000 极值点：(x,y) = (1.000000,1.000000) 极值：f(x,y) = 0.00000000000000000000   x_best =        1      1     f_best =   0.0
　　实例4
　　主程序 clc; clear all; close all; syms x; f =  9.*x.^2-sin(x)-1; [x_optimization,y] = Newton_Method(f,2); x_optimization = double(x_optimization); y =vpa(y) x_optimization x = -10:0.01:10; ft = 9.*x.^2-sin(x)-1; figure(1) plot(x,ft); hold on; plot(x_optimization,y,＂r*＂);
　　Newton_Method函数程序 function [x_optimization,f_optimization] = Newton_Method(f,x0) format long; %   f：目标函数 %   x0：初始点 %   epsilon：精度 %   x_optimization：目标函数取最小值时的自变量值 %   f_optimization：目标函数的最小值 if nargin == 2     epsilon = 1.0e-6; end df = diff(f);       %   一阶导数 d2f = diff(df);     %   二阶导数 k = 0; dfxk = 1; xk = x0; while dfxk > epsilon     dfx = subs(df,symvar(df),xk);     if diff(d2f) == 0         d2fx = double(d2f);     % 二阶导数不能为零     else         d2fx = subs(d2f,symvar(d2f),xk);      end     xk_next = xk - dfx/d2fx;         k = k + 1;                       dfxk = abs(dfx);     xk = xk_next;   %   迭代 end   x_optimization = xk_next; f_optimization = subs(f,symvar(f),x_optimization); format short; end
　　运行结果   y =   -1.0277492701423876507411151284973     x_optimization =       0.0555
　　本文内容来源于网络，仅供参考学习，如内容、图片有任何版权问题，请联系处理，24小时内删除。
　　作 者 | 郭志龙
　　编 辑 | 郭志龙
　　校 对 | 郭志龙