今天主要是讲解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 yf = 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_bestf_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;endx0 = 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\n",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.65572826847430654151x_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 yfx = 4*(x-y)-x^2-y^2;%定义二元变量函数x0 = [1 1];%初始点 x0(1,1)[x_best,f_best] = Newton(fx,x0,[x y]);x_bestf_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.00000000000000000000x_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 yf = (1-x)^2+100*(y-x^2)^2;x0 = [0 0];%初始点 x0(1,1)[x_best,f_best] = Newton(f,x0,[x y]);x_bestf_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.00000000000000000000x_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_optimizationx = -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;enddf = 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; % 迭代endx_optimization = xk_next;f_optimization = subs(f,symvar(f),x_optimization);format short;end运行结果
y = -1.0277492701423876507411151284973 x_optimization = 0.0555本文内容来源于网络，仅供参考学习，如内容、图片有任何版权问题，请联系处理，24小时内删除。
作 者 | 郭志龙
编 辑 | 郭志龙
校 对 | 郭志龙
