最速下降法MATLAB程序实现.docx

function [drag,ddrag]=fxy(a0,s0)%要优化的目标函数%返回值：drag，函数在（a0，s0）点处的值；ddrag，函数在（a0，s0，drag）处的梯度（行向量）%参数：（a0,s0） function d = f(a,s) w = 0.12754*a^1.5 + ( 0.12754^2*a^3 + 4*(4940/s + 45.42) )^0.5; w = w^2*s/4; cl = 2*w / (1506.75*s); cd = 0s + 0* a^0.1 * s^-0.1 + cl^2*a^-1/(pi*0.96); d = w*cd/cl; endddrag = zeros(1,2);drag = f(a0,s0);h = 1e-6;ddrag(1) = (f(a0+h,s0)-drag)/h; %函数对a的偏导数ddrag(2) = (f(a0,s0+h)-drag)/h; %函数对s的偏导数end%最速下降法%-------------------------------------------------%-------------------------------------------------clear; %绘制等高线x=5:0.1:30; %展弦比Ay=5:0.1:30; %机翼面积S[X,Y]=meshgrid(x,y);dim = size(X);for i=1:dim(1) for j=1:dim(2) [Z(i,j),null]=fxy(X(i,j),Y(i,j)); endend%mesh(X,Y,Z);hold on;[c,h]=contourf(X,Y,Z,10);clabel(c,h);colormap coolxlabel(A);ylabel(S);%-------------------------------------------------%-------------------------------------------------eg = 1e-6; %精度要求ea = 1e-6;er = 1e-2;num = 0; %迭代次数num1 = 0; %确定最优步长迭代的次数x0=15;y0=30; %起始点[z0,dz0] = fxy(x0,y0);while norm(dz0)eg z00 = z0; x00 = x0; y00 = y0; h = norm(dz0)/1e4; %初始步长 dz0 = dz0/norm(dz0); %沿梯度方向的单位向量 x1 = x0 - dz0(1)*h; y1 = y0 - dz0(2)*h; [z1,dz1] = fxy(x1,y1); step0 = 0; %[step0,step1]，步长极小值区间 step1 = step0+h; step = step0; if z1z0 while z1z0 %用进退法确定步长的极值区间 step = step0; step0 = step1; z0 = z1; h = 2*h; step1 = step1 + h; x1 = x0 - dz0(1)*step1; y1 = y0 - dz0(2)*step1; [z1,dz1] = fxy(x1,y1); end step1 = step0; step0 = step;end while (step1-step0)eg %黄金分割法求最优步长 x2 = x0 - dz0(1)*( step0 + 0.618*(step1-step0) ); y2 = y0 - dz0(2)*( step0 + 0.618*(step1-step0) ); [z2,dz2] = fxy(x2,y2); x1 = x0 - dz0(1)*( step0 + 0.382*(step1-step0) ); y1 = y0 - dz0(2)*( step0 + 0.382*(step1-step0) ); [z1,dz1] = fxy(