makispaiktis

Texnikes_Ergasia2_Thema4_ZitoumenoC

Dec 5th, 2020
952
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
  1. clear all
  2. clc
  3.  
  4. xxx = -1;
  5. yyy = -1;
  6. elaxistoMeGammaMetavlhto(xxx, yyy, 0.8, 0);
  7. elaxistoMeGammaMetavlhto(xxx, yyy, 1, 1);
  8.  
  9.  
  10. function elaxistoMeGammaMetavlhto(x0, y0, s, flag)
  11.    
  12.     % 1. Ορίζω το ε της συνθήκης τερματισμού ίσο με 2/1000
  13.     epsilon = 0.002;
  14.     syms x y
  15.     f = x^3 * exp(-x^2-y^4);
  16.     klisi = gradient(f, [x,y]);
  17.     KLISI = matlabFunction(klisi);
  18.     essianos = jacobian(klisi)
  19.     ESSIANOS = matlabFunction(essianos);        % Δίνει έναν 2x2 πίνακα
  20.    
  21.     % 2. Ορίζω τις λίστες που θα τοποθετήσω τα xi, yi. Τις ονομάζω xList,
  22.     % yList, θα βάζω επίσης και τις τιμές της f (fList) και του μέτρου της
  23.     % κλίσης της (normKlisisList) σε άλλες 2 λίστες
  24.     k = 1;
  25.     xList = [];             yList = [];
  26.     xList(1) = x0;          yList(1) = y0;
  27.     fList = [];             normKlisisList = [];        gammaList = [];
  28.     fList(1) = subs(f, {x,y}, {xList(length(xList)), yList(length(yList))});
  29.     normKlisisList(1) = norm(subs(klisi, {x,y}, {xList(length(xList)), yList(length(yList))}));
  30.    
  31.     % Armijo Constants
  32.     a = 1/1000;             b = 0.3;        
  33.     gammaList(1) = s;
  34.    
  35.     while normKlisisList(length(normKlisisList)) > epsilon
  36.         k = k + 1;
  37.         % *********************************************************************************
  38.         % ************************ Levenberg - Marquardt **********************************
  39.         % *********************************************************************************
  40.        
  41.         % Πρέπει να βρω αν ο εσσιανός στο xk είναι θετικά ορισμένος, θα
  42.         % γίνει μέσω ιδιοτιμών
  43.         ESS = ESSIANOS(xList(k-1), yList(k-1))
  44.         EIG = eig(ESS);
  45.         MIN = min(EIG);
  46.         dk = [0;0];                 % A 2x1 vector
  47.        
  48.         if MIN > 0
  49.             % Τότε, η ελάχιστη ιδιοτιμή του εσσιανού στο σημείο xk είναι
  50.             % θετική και άρα ο εσσιανός είναι θετικά ορισμένος
  51.             disp('Thetika orismenos');
  52.             inv_essianos = inv(ESS)
  53.             dk = - inv_essianos * KLISI(xList(k-1), yList(k-1));        % A 2x1 vector
  54.         else
  55.             % Αν η ελάχιστη ιδιοτιμή του είναι -4 π.χ. πρέπει να προσθε΄σω
  56.             % 4 + (κάτι) και όλο αυτό επί τον μοναδιαίο, π.χ. 4,4Ι
  57.             disp('Oxi thetika orismenos');
  58.             mk = - MIN * 1.1
  59.             A = ESS + mk * eye(2);
  60.             inv_A = inv(A)
  61.             % Τώρα, ο Α είναι σίγουρα θετικά ορισμένος και αντιστρέφεται
  62.             dk = - inv_A * KLISI(xList(k-1), yList(k-1));
  63.         end
  64.        
  65.         % ********************************************************************************************
  66.         % ******************************* Κανόνας Armijo *********************************************
  67.         % ********************************************************************************************
  68.         % Δοκιμάζω τιμές ακέραιες 0,1,2,... στο mk και φτιάχνω το βήμα γk
  69.         % και κατ' επέκταση και το xk+1
  70.         for mk = 0:10
  71.             k
  72.             mk
  73.             gk = s * b^mk;
  74.             xNew = xList(k-1) + gk * dk(1);
  75.             yNew = yList(k-1) + gk * dk(2);
  76.             F = matlabFunction(f);
  77.             DF = matlabFunction(klisi);
  78.             % Ήρθε η ώρα να τσεκάρουμε την ανίσωση
  79.             aristeroMelos = F(xList(k-1), yList(k-1)) - F(xNew, yNew);
  80.             grad = DF(xList(k-1), yList(k-1));
  81.             dexioMelos = a * b^mk * s * grad' * grad;
  82.             if aristeroMelos > dexioMelos
  83.                 gammaList(k) = gk;
  84.                 xList(k) = xNew;
  85.                 yList(k) = yNew;
  86.                 fList(k) = F(xNew, yNew);
  87.                 normKlisisList(k) = norm(subs(klisi, {x,y}, {xList(length(xList)), yList(length(yList))}));
  88.                 display('Success for this step k with mk in value of:')
  89.                 mk
  90.                 break;
  91.             end
  92.         end
  93.            
  94.            
  95.     end
  96.    
  97.     xList
  98.     yList
  99.     fList
  100.     normKlisisList
  101.     k
  102.     gammaList
  103.     % Τα τελευταία xk, yk κάθε λίστας τα ονομάζω εν συντομία xx και yy
  104.     display('~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~')
  105.     display('~~~~~~~~ About the last found xk (xx) and yk (yy) ~~~~~~~~')
  106.     display('~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~')
  107.     xx = xList(length(xList))
  108.     yy = yList(length(yList))
  109.     F_xx_yy = fList(length(fList))
  110.     NORM_KLISIS = normKlisisList(length(normKlisisList))
  111.     if flag == 0
  112.         plot(fList, 'ro');
  113.         title('Red for s = 0.8, blue for s = 1');
  114.         xlabel('Steps k');
  115.         ylabel('Values of f');
  116.         hold on
  117.     else
  118.         plot(fList, 'bx');
  119.         title('Red for s = 0.8, blue for s = 1');
  120.         xlabel('Steps k');
  121.         ylabel('Values of f');
  122.         hold on
  123.     end
  124.     display('**********************************************************')
  125. end
RAW Paste Data