本文介绍了几种数值优化方法，包括二分法、共轭梯度法、内点罚函数法、最速下降法和牛顿法，并提供了相应的MATLAB代码示例。此外，还展示了绘制统计图的过程，特别是柱状图与折线图的组合展示。
文章目录
一、二分法
用二分法求函数y=x^2-sin(x) 的极小值，绘制出迭代过程，自己选定初始迭代区间，要求最终精度优于0.01。
format short
s=input('请输入函数表达式：f = ','s');
f=inline(s);
a = input('请输入区间左端点:a=');
b = input('请输入区间右端点:b=');
eps=input('请输入停止精度要求:eps='); %%“|b-x|<=eps”
k=1;
x=(a+b)/2;
fprintf(' k a f(a) b f(b) xk f(xk)\n ');
T=[k,a,f(a),b,f(b),x,f(x)];
while abs(T(k,4)-T(k,6))>eps/2
k=k+1;
if f(x)*f(a)==0
a=a;
b=x;
x=(a+b)/2;
T=[T;k,a,f(a),b,f(b),x,f(x)];
break
elseif f(x)*f(a)>0
a=x;
b=b;
x=(a+b)/2;
T=[T;k,a,f(a),b,f(b),x,f(x)];
elseif f(x)*f(a)<0
a=a;
b=x;
x=(a+b)/2;
T=[T;k,a,f(a),b,f(b),x,f(x)];
end
end
disp(T);
fprintf('经过%d次迭代，函数方程根的近似解为：x=%.8f\n',k-1,T(k-1,6))%利用二分法，求解函数f(x)=x×sinx -x/6在[0,1]上的最小值，要求精度小于0. 1，给出每一步的详细计算结果；
a0=0;b0=1;i=0;%初始区间
x=0;
while (b0-a0)>0.1
x=(a0+b0)/2;i=i+1;
f = sin(x) + x*cos(x)-1/6;
if f>0%导数大于零，令区间右侧为x
b0=x;
disp("第"+i+"区间为"+"["+a0+","+b0+"]");
end
if f<0%导数小于零，令区间左侧为x
a0=x;
disp("第"+i+"区间为"+"["+a0+","+b0+"]");
end
if f==0%导数等于零，结束
break;
end
end
disp("最终结果为"+"["+a0+","+b0+"]"+"一共迭代了"+i+"次"); 二、共轭梯度法
代码如下（示例）：
syms x1 x2 x3;
xk=[0;0;0];
Q=[3 0 1;0 4 2;1 2 3];
x=[x1;x2;x3];
b=[3;0;1];
f=1/2*x'*Q*x-b'*x;
grad_f=jacobian(f,x);
grad_f=transpose(grad_f);
flag=double(subs(grad_f,x,xk));
dk=-flag;
z=zeros(3,1);
i=0;
while i<=3
if isequal（flag，z）
a=-(dk'*flag)/(transpose(dk)*Q*dk);
a=double(a);
xk=xk+a*dk;
flag=subs(grad_f,x,xk);
flag=double(flag);
i=i+1;
if isequal（flag，z）
b=(dk'*Q*flag)/(dk'*Q*dk);
b=double(b);
dk=-flag+b*dk;
dk=double(dk);
else
x_final=xk;
disp(x_fianl);
break;
end
end
end
三、内点罚函数法
代码如下（示例）：
用内点罚函数法求解问题，内点罚函数为
x0=[4 4]';
x1 = linspace(-4,4);
y1 = linspace(-4,4);
[X,Y] = meshgrid(x1,y1);
Z = X.^4 + 2*Y.^4;
figure(1),[C,h]=contour3(X,Y,Z);
clabel(C,h);hold on;
maxk=5000; %最大迭代次数
rho=0.5;
sigma=0.4;
k=0;
e=0.1; %精度
while(k<maxk)
g=gfun(x0); %计算梯度
d=-g;
if(norm(d)<e),break;end
%用Amrijo搜索技术确定步长
m=0;mk=0;
while(m<20) %最大迭代次数
if(fun(x0+rho^m*d)<fun(x0)+sigma*rho^m*g'*d)
mk=m;
break;
else
m=m+1;
end
end
xx=x0;
x0=x0+d*rho^mk;
fprintf('第%d次迭代：x =%f y=%.3f',k+1,x0(1),x0(2));
% disp(x0(1),x0(2));
vall=fun(xx);
val=fun(x0);
disp([' 此时: f(x) = ',num2str(val)])
figure(1),plot3(x0(1),x0(2),val, 'r.','markersize',10);
figure(1),line([xx(1),x0(1)],[xx(2),x0(2)],[vall,val]);
k=k+1;
end
x=x0;
val=fun(x0);
fprintf('最优解为：x =%f y=%.3f',x0(1),x0(2));
disp([' 此时: f(x) = ',num2str(val)])
四、最速下降法
代码如下（示例）：
绘制
的四幅等高线图作为子图，并将四幅子图排成2行2列
% 定义变量取值范围
x1 = linspace(-1, -0, 500)
x2 = linspace(0, 1, 500);
% 创建网格矩阵
[X, Y] = meshgrid(x1, x2);
% 定义内点罚函数
penalty_function = @(x1, x2, u) x1.^2 + x2.^2 - u * log(-x1 + x2 - 1);
% 定义内点罚函数参数
u = 0.1
% 计算内点罚函数值
Z = penalty_function(X, Y, u);
% mesh(X,Y,real(Z));
% 绘制等高线图
subplot(221);
contour(X, Y, real(Z), 20);
title('u=0.100');
xlabel('x1');
ylabel('x2');
subplot(222);
u=0.050;
Z = penalty_function(X, Y, u);
% surf(X,Y,real(Z));
contour(X, Y, real(Z), 20);
title('u=0.050');
xlabel('x1');
ylabel('x2');
subplot(223);
u=0.010;
Z = penalty_function(X, Y, u);
contour(X, Y, real(Z), 20);
title('u=0.010');
xlabel('x1');
ylabel('x2');
subplot(224);
u=0.001;
Z = penalty_function(X, Y, u);
contour(X, Y, real(Z), 20);
title('u=0.001');
xlabel('x1');
ylabel('x2');五、牛顿法
代码如下（示例）：
%利用Newton法，求解函数f(x)=x×sinx -x/6在[0,1]上的最小值，要求精度小于0. 1，给出每一步的详细计算结果；
a0=0;b0=1;i=0;%初始区间
x=0.5;t=1;
while abs(t-x)>0.1
f1=sin(x)+x*cos(x)-1/6;%一阶导
f2=cos(x)+cos(x)-x*sin(x);%二阶导
t=x;%区间值修改
x=x-(f1/f2);%x值更改
i=i+1;%迭代次数加1
disp("第"+i+"次迭代结果为："+x+",f（"+x+"）="+f);
end
disp("最终结果为:x="+x+",y="+f+"一共迭代了"+i+"次");
六、割线法
eps=5e-8;%割线法
delta=1e-10;
N=100;%迭代次数上限
k=0;x0=0;x1=1;%给定两个点，分部在根的两侧
tic
while(1)
x2=x1-func2_2(x1)*(x1-x0)/(func2_2(x1)-func2_2(x0));
k=k+1
if(k>N||abs(x2)<eps)
disp('Newton method failed')
break
end
if abs(x2)<1
d=x2-x1;
else
d=(x2-x1)/x2;
end
x0=x1;
x1=x2
y1=func2_2(x1)
if(abs(d)<eps||abs(func2_2(x2))<delta)
break
end
end
toc
function y=func2_2(x)%函数
y=2*exp(-x)-sin(x);
end七、绘制统计图
%定义线上点的x坐标
x = 1:1:4 %第一个1是起始端点，中间的是步长，最后一个是结束端点
%纵坐标
samp1 = [99.00 95.00 94.00 94.00]; %A
samp3 = [30.00 40.00 50.00 60.00]; %C
fig = figure('Units','centimeter','Position',[5 5 12 13]);
%中括号里面的是调整图片大小的，主要调整最后面两个长度和宽度
left_color = [153/255,70/255,180/255];%左边y坐标的颜色
right_color = [153/255,154/255,205/255]; %右边y坐标的颜色
set(fig,'defaultAxesColorOrder',[left_color; right_color]);
hold on;
% 画柱状图
yyaxis left
bar1(:,1) = samp1;
GO = bar(bar1,0.5,'EdgeColor','black');%边框颜色为黑色
GO.FaceColor = 'flat';%设置柱状图的颜色
GO.CData(1,:)=[140/255,94/255,141/255];
GO.CData(2,:)=[228/255,191/255,203/255];
GO.CData(3,:)=[204/255,129/255,152/255];
GO.CData(4,:)=[156/255,136/255,179/255];
ylim([0 110]) %y轴显示的范围，根据需要调整
ylabel('Price(yuan)','FontName', 'Times New Roman')
% 标记数据到柱状图
offset_vertical = 2.5; % 高度
offset_horizon = 0; % 左右偏离
for i = 1:length(samp1)
if samp1(i)>=0
text(i - offset_horizon,samp1(i) + offset_vertical,num2str(samp1(i),'%.2f'),'FontName', 'Times New Roman','VerticalAlignment','middle','HorizontalAlignment','center');
else
text(i - offset_horizon,samp1(i) - offset_vertical,num2str(samp1(i),'%.2f'),'FontName', 'Times New Roman','VerticalAlignment','middle','HorizontalAlignment','center');
end
end
% 右边y轴（折线图）
yyaxis right
Line1=plot(x,samp3,'*--','LineWidth', 1.5,'MarkerEdgeColor', '#c2d7f3');
%MarkerEdgeColor 标记物颜色;
Line1.Color=[110/255, 123/255, 178/255];
ylim([0 110]) %右边y轴显示范围
ylabel('Cost(yuan)','FontName', 'Times New Roman')
% 标记数据到折线图
offset_vertical = 2.5; % 高度
offset_horizon = 0; % 左右偏离
for i = 1:length(samp3)
if samp3(i)>=0
text(i - offset_horizon,samp3(i) + offset_vertical,num2str(samp3(i),'%.2f'),'FontName', 'Times New Roman','VerticalAlignment','middle','HorizontalAlignment','center');
else
text(i - offset_horizon,samp3(i) - offset_vertical,num2str(samp3(i),'%.2f'),'FontName', 'Times New Roman','VerticalAlignment','middle','HorizontalAlignment','center');
end
end
% 图例
legend({'Price(yuan','Cost(yuan)'});
title("Commodity");
% x轴
grid on; %画网格
xticks([1 2 3 4])
xticklabels({'commodity1','commodity2','commodity3','commodity4'})
扫一扫
2万+
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
