makispaiktis

Texnikes_Ergasia2_Thema3_ZitoumenoC

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