最速下降法Matlab程序

%最速下降梯度法matlab程序 % Steepest Descent Method % By Kshitij Deshpande clc clear all warning off prompt = {\'Coeficients if X1=\',\'Coefficients of X2=\',\'Coefficeint of X1X2=\',\'Initial Point=\'}; def = {\'[2 1 0]\',\'[1 -1 0]\',\'2\',\'[0 0]\'}; a=inputdlg(prompt,\'Data\',1,def); a=char(a); [m,n]=size(a); x1 = eval(a(1,1:n));x2=eval(a(2,1:n));x1x2=eval(a(3,1:n));X1=eval(a(4,1:n)); delf1(1) = polyval(polyder(x1),X1(1)); delf1(1) = (delf1(1))+(x1x2*X1(2)); delf1(2) = polyval(polyder(x2),X1(1)); delf1(2) = (delf1(2))+(x1x2*X1(1)); s=-delf1; %%%%%%%%%% %report srep(1,1:2)=s; %%%%%%%%%% x1new(1)=s(1)^2;x1new(2)=2*X1(1)*s(1);x1new(3) = X1(1)^2; x1new=x1new*x1(1); x1new_(2)=x1(2)*s(1);x1new_(3)=x1(2)*X1(1); x1new = x1new+x1new_; x2new(1)=s(2)^2;x2new(2)=2*X1(2)*s(2);x2new(3) = X1(2)^2; x2new=x2new*x2(1); x2new_(2)=x2(2)*s(2);x2new_(3)=x2(2)*X1(2); x2new = x2new+x2new_; x1x2new(1)=s(1)*s(2);x1x2new(2)=X1(1)*s(2)+X1(2)*s(1);x1x2new(3)=X1(1)*X1(2); x1x2new=x1x2*x1x2new; df = polyder(x1new+x2new+x1x2new); lambda(1) = roots(df); X1=X1+lambda(1)*s; Xrep(1,1:2)=X1; delf1(1) = polyval(polyder(x1),X1(1)); delf1(1) = (delf1(1))+(x1x2*X1(2)); delf1(2) = polyval(polyder(x2),X1(2)); delf1(2) = (delf1(2))+(x1x2*X1(1)); if all(X1)== 0 fprintf(\'%d %d is the optimum point\',X1(1),X1(2)); end itrep(1)=1; it=2; while all(delf1)==1 s=-delf1; x1new(1)=s(1)^2;x1new(2)=2*X1(1)*s(1);x1new(3) = X1(1)^2; x1new=x1new*x1(1); x1new_(2)=x1(2)*s(1);x1new_(3)=x1(2)*X1(1); x1new = x1new+x1new_; x2new(1)=s(2)^2;x2new(2)=2*X1(2)*s(2);x2new(3) = X1(2)^2; x2new=x2new*x2(1); x2new_(2)=x2(2)*s(2);x2new_(3)=x2(2)*X1(2); x2new = x2new+x2new_; x1x2new(1)=s(1)*s(2);x1x2new(2)=X1(1)*s(2)+X1(2)*s(1);x1x2new(3)=X1(1)*X1(2); x1x2new=x1x2*x1x2new; df = polyder(x1new+x2new+x1x2new); lambda(it) = roots(df); X1=X1+lambda(it)*s; delf1(1) = polyval(polyder(x1),X1(1)); delf1(1) = (delf1(1))+(x1x2*X1(2)); delf1(2) = polyval(polyder(x2),X1(2)); delf1(2) = (delf1(2))+(x1x2*X1(1)); itrep(it)=it; srep(it,1:2)=s; Xrep(it,1:2)=X1; it=it+1; end [m,n]=size(itrep); matrix=[itrep\' srep(1:n,1) srep(1:n,2) Xrep(1:n,1) Xrep(1:n,2)]; answer = char(num2str(X1)); answer = [\'The optimal point is [\' answer \']\']; msgbox(answer,\'Solution\'); disp(\' Press Any key to View Detailed Report............\'); pause echo off report steep; clc -------------------------------------------------------------------------------- %最速下降法（爬山法）的一个matlab程序 function y=steepest(x) %This program uses the steepest descent direction algorithm %to calculate the minimum of the function f(x)=x(1)^2+2*x(2)^2 format long eps=input(\'please input your accuracy:\'); %eps is the demmanded accuracy on the norm of %the gradient of the objective function m=1; %m is the count of the iteration step of the algorithm iterstep(1,:)=x; %iterstep contains the intermediate points of iteration while norm(gradobject1(x))>eps grad=gradobject1(x); alpha=goldsplictobj(x); x=x-alpha*grad; iterstep(m+1,:)=x; m=m+1; end step=max(size(iterstep))-1 plot(iterstep(:,1),iterstep(:,2)); %Draw the search trajectory title(\'The search trajectory of the Steepest dscent direction algorithm\'); xlabel(\'x1-axis\'); ylabel(\'x2-axis\'); text(x(1),x(2),\'The minimum point found by the algorithm\'); text(iterstep(1,1),iterstep(1,2),\'The initial point (2,1)\'); gtext(\'The number of the total iteration steps of the algorithm is:\'); gtext(\'The set accuracy in advance is 1.0*10^{-10}\'); %The following subfunction is on the objective function function y=object1(v) y=v(1)^2+2*v(2)^2; %The following subfunction is on the gradient of %the objective function function y=gradobject1(v) y(1)=2*v(1); y(2)=4*v(2); %The following subfunction is on the comming %search function of alpha function y=substi(alpha,x) y=feval(\'object1\',x-alpha*gradobject1(x)); %The following subfunction is on the goldspliction %search of the substi function function y=goldsplictobj(x) a=0; b=10; eps=0.01; y1=a+0.382*(b-a); y2=a+0.618*(b-a); while abs(b-a)>eps if substi(y1,x)>substi(y2,x) a=y1; b=b; y1=a+0.382*(b-a); y2=a+0.618*(b-a); elseif substi(y2,x)>substi(y1,x) a=a; b=y2; y1=a+0.382*(b-a); y2=a+0.618*(b-a); else a=y1; b=y2; y1=a+0.382*(b-a); y2=a+0.618*(b-a); end end y=(y1+y2)/2;