Steepest Descent

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Uebungsaufgabe Quadratische Funktionen mit dem Gradientenverfahren
% Copyright 2007 Thomas Ludwig
% This file is distributed under the terms of the GNU General Public License,
% see http://www.gnu.org/copyleft/gpl.txt
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Mit dieser Hauptfunktion wird die Berechnung gestartet.
% Der Parameter c ruft vorbereitete Beispiele in der Funktion
% gradient_beispiele auf.
function gradienten( c );

    % Uebung Gradientenverfahren
    % Uebungsaufgabe: quadratische Funktionen
    % f(x,y) = a x^2 + b x y + c y^2 + d x + e y + f
    % Bestimmen Sie das Minimum von f und dessen Lage
   
    % Startwerte
    x0=[0, 0];
   
   % Beispiele auswaehlen
   fx= gradient_beispiele(c);

   % Der Numerische Lösungsansatz
   func_rechenknecht= @(x) rechenknecht(x, fx);
   [fx.xy,fx.fval,fx.exitflag,fx.output,fx.grad,fx.hessian] = fminunc(func_rechenknecht, x0, fx.options);
  
   % Der Grafische Lösungsansatz
   grafik(fx);

   % Ausgabe
    disp('Numerische Lösung:');
    fprintf('Lokales Extrema x = [%g %g] \n\n',fx.xy(1),fx.xy(2))
    disp(fx.output)  
end

% Rechenknecht beinhaltet die Funktion und den Gradienten
function [f,g]= rechenknecht(x, fx);        % anonymous functions
    f= fx.a*x(1).^2 + fx.b*x(1)*x(2) + fx.c*x(2).^2 + fx.d*x(1) + fx.e*x(2) + fx.f;     % Quadratische Funktion
    g= [2*fx.a*x(1) + fx.b*x(2) + fx.d; 2*fx.c*x(2) + fx.b*x(1) + fx.e];             % Gradient
    H= [2*fx.a, fx.b; fx.b, 2*fx.c];                                                 % Hesse/Hessian
   
    % Muss Maximiert werden?
    if ( min(eig(H)) < 0 );
        g= -g;
        disp('Diese Funktion muss Maximiert werden -> siehe Grafische Lösung');
    end 
end

%
% Grafische Lösung
%
function grafik(fx);
  
    schrittweite= 0.25;
    xk= -2:schrittweite:2;
    yk= -2:schrittweite:2;
    [Xk Yk] = meshgrid(xk,yk);

    Zk= fx.a*Xk.^2 + fx.b*Xk*Yk + fx.c*Yk.^2 + fx.d*Xk + fx.e*Yk;
   
    % Bestimmt den Gradient der Funktion Z
    [Fx,Fy]= gradient(Zk,schrittweite,schrittweite);

    % Zeichnet das Richtungsfeld des Gradienten
    subplot(3,2,1:2);
        % Erzeugt die Hoehenlinien zu x, y und der Funktion z
        contour(xk,yk,Zk), hold on;
       
        % Zeichnet die Richtung des Gradienten ein
        quiver(Xk,Yk,Fx,Fy), hold off;
        xlabel('x'), ylabel('y'), title('Grafische Lösung');

    subplot(3,2,3);
        contourf(xk,yk,Zk), colorbar('horiz'), hold on;
        quiver(Xk,Yk,Fx,Fy), hold off;
        xlabel('x'), ylabel('y'), zlabel('z');
       
    subplot(3,2,4);  
        mesh(xk,yk,Zk);
        xlabel('x'), ylabel('y')
       
    subplot(3,2,5);  
        surfc(Xk,Yk,Zk);
        xlabel('x'), ylabel('y');
       
    subplot(3,2,6);
        % Normierung des Gradienten auf Laenge 1 und Zeichnen
        Fxn=max(abs(Fx),0.0000000001);  % Damit nicht durch Null geteilt wird,
        Fyn=max(abs(Fy),0.0000000001);  % an Nullstellen des Gradienten.
        Fn=sqrt(Fxn.^2+Fyn.^2);
        quiver(Xk,Yk,Fx./Fn,Fy./Fn);
        xlabel('x'), ylabel('y'), grid on;
end

%
% Koeffizienten fuer verschiedene Berechnungen
% Parameter c waehlt die gewuenschte Aufgabe aus
%
function fx= gradient_beispiele(c);

    switch c
       
        case 1
            % x² und y² sind hier positiv
            % Somit ist die Funktion nach oben geoeffnet und besitzt
            % ein Minimum. Deshalb suchen wir ein lokales Minimum    
            % Dies ist dann eine Quelle (das lokale Minimum)
            fx.a= 1.5;
            fx.b= 2;
            fx.c= 1;
            fx.d= -1;
            fx.e= -2; 
            fx.f= 0;
            fx.Z= inline('1.5*x(1).^2 + 2*x(1)*x(2) + 1*x(2).^2 -1*x(1) - 2*x(2)');
           
        case 2
            % x² und y² sind hier negativ
            % Somit ist die Funktion nach unten geoeffnet und besitzt
            % ein Maximum. Deshalb suchen wir ein lokales Maximum
            fx.a= -2.5;
            fx.b= -4;
            fx.c= -3;
            fx.d= -4;
            fx.e= -2.5; 
            fx.f= 0;
            fx.Z= inline('-2.5*x(1).^2 - 4*x(1)*x(2) - 3*x(2).^2 - 4*x(1) - 2.5*x(2)');
           
        case 3
            % x² und y² sind hier positiv
            % Somit ist die Funktion nach oben geoeffnet und besitzt
            % ein Minimum. Deshalb suchen wir ein lokales Minimum    
            % Dies ist dann eine Quelle (das lokale Minimum)
            fx.a= 1.5;
            fx.b= 2;
            fx.c= 1;
            fx.d= 0.5;
            fx.e= 1; 
            fx.f= 0;
            fx.Z= inline('1.5*x(1).^2 - 2*x(1)*x(2) - 1*x(2).^2 - 0.5*x(1) - 1*x(2)');
           
        case 4
            % x² und y² sind hier negativ
            % Somit ist die Funktion nach unten geoeffnet und besitzt
            % ein Maximum. Deshalb suchen wir ein lokales Maximum
            fx.a= -2.5;
            fx.b= 4;
            fx.c= -3;
            fx.d= 2;
            fx.e= 1.25; 
            fx.f= 0;     
            fx.Z= inline('-2.5*x(1).^2 + 4*x(1)*x(2) - 3*x(2).^2 + 2*x(1) + 1.25*x(2)');
    end

    % Optionen fuer fminunc
    clear options;
    options= optimset('LargeScale','off', 'Display','off');  % LargeScale off wegen Hesse
    options= optimset(options,'HessUpdate','steepdesc');     % Hesse mit steepest descent benutzen
    options= optimset(options,'LineSearchType','quadcubic'); % Quadratic polynomial interpolation and extrapolation method
    options= optimset(options,'gradobj','on');               % Gradienten Objekt benutzten
    fx.options= optimset(options, 'TolFun',1e-10, 'TolX',1e-10, 'MaxIter',1000, 'MaxFunEvals',250); 
  
end