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小时内删除。

作 者 | 郭志龙

编 辑 | 郭志龙
校 对 | 郭志龙