makispaiktis

Texnikes_Ergasia2_Thema3_ZitoumenoB

Dec 1st, 2020
660
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.1, 1, 0);
  5. elaxistoMeGammaMetavlhto(-1, -1, 0.5, 1.5, 1);
  6.  
  7.  
  8. function elaxistoMeGammaMetavlhto(x0, y0, akro1, akro2, 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.     gammaList(1) = 0;
  29.    
  30.     while normKlisisList(length(normKlisisList)) > epsilon
  31.         k = k + 1;
  32.         % dk = -subs(klisi, {x,y}, {xList(length(xList)), yList(length(yList))});     % η κλίση πριν
  33.         % Το dk είναι 2*1 διάνυσμα, το 1ο στοιχείο του αφορά τον υπολογισμό
  34.         % του x και το 2ο τον υπολογισμό του y
  35.         % Πρέπει το γ που θα επιλέξω να κάνει minimize την f(xk + γkdk)
  36.         xPrin = xList(k-1);
  37.         yPrin = yList(k-1);
  38.         dkVector = DK(xList(k-1), yList(k-1));          % A 2x1 vector
  39.        
  40.         % ********************************************************************************************
  41.         % ************************** Εσωτερική Βελτιστοποίηση ****************************************
  42.         % ********************************************************************************************
  43.         gamma = internalOptimization(f, xPrin, yPrin, dkVector, akro1, akro2);
  44.        
  45.        
  46.         xList(k) = xPrin + gamma * dkVector(1);
  47.         yList(k) = yPrin + gamma * dkVector(2);
  48.         fList(k) = subs(f, {x,y}, {xList(length(xList)), yList(length(yList))});
  49.         normKlisisList(k) = norm(subs(klisi, {x,y}, {xList(length(xList)), yList(length(yList))}));
  50.         gammaList(k) = gamma;
  51.     end
  52.    
  53.     xList
  54.     yList
  55.     fList
  56.     normKlisisList
  57.     k
  58.     gammaList
  59.     % Τα τελευταία xk, yk κάθε λίστας τα ονομάζω εν συντομία xx και yy
  60.     disp('~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~')
  61.     display('~~~~~~~~ About the last found xk (xx) and yk (yy) ~~~~~~~~')
  62.     display('~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~')
  63.     xx = xList(length(xList))
  64.     yy = yList(length(yList))
  65.     F_xx_yy = fList(length(fList))
  66.     NORM_KLISIS = normKlisisList(length(normKlisisList))
  67.     if flag == 0
  68.         plot(fList, 'ro');
  69.         title('Red for gamma = 0.2, blue for gamma = 0.5');
  70.         xlabel('Steps k');
  71.         ylabel('Values of f');
  72.         hold on
  73.     else
  74.         plot(fList, 'bx');
  75.         title('Red for 0.1 <= gamma <= 1, blue for 0.5 <= gamma <= 1.5');
  76.         xlabel('Steps k');
  77.         ylabel('Values of f');
  78.         hold on
  79.     end
  80.     display('**********************************************************')
  81. end
  82.  
  83.  
  84.  
  85.  
  86.  
  87.  
  88.  
  89.  
  90. function gamma = internalOptimization(f, xPrin, yPrin, dk, akro1, akro2)
  91.     syms x y G
  92.     X = xPrin + G * dk(1);
  93.     Y = yPrin + G * dk(2);
  94.     % Πλέον, τα νέα X,Y έχουν σαν μόνη άγνωστη μεταβλητή το G = γk
  95.     % Η F που έχει προκύψει είναι μόνο συνάρτηση του G
  96.     F = subs(f, {x,y}, {X, Y});
  97.     DF = diff(F, 'G');
  98.     l = 0.005;
  99.     counter = 0;
  100.     while akro2 - akro1 > l
  101.         counter = counter + 1;
  102.         kentro = (akro1 + akro2) / 2;
  103.         paragwgos = subs(DF, kentro);
  104.         if paragwgos > 0
  105.             akro2 = kentro;
  106.         else
  107.             akro1 = kentro;
  108.         end
  109.      end
  110.      % Οι μεταβλητές akra1, akra2 μετά το τέλος της διαδικασίας έχουν
  111.      % κάποιες διαμορφωμένες τιμές. Θα κάνω τον μέσο όρο τους για το τελικό
  112.      % γ που θα επιλέξω να επιστρέψω
  113.      gamma = (akro2 + akro1) / 2;
  114. end
RAW Paste Data