今天主要是讲解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\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.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小时内删除。
作 者 | 郭志龙
编 辑 | 郭志龙
校 对 | 郭志龙