makispaiktis

Texnikes_Ergasia2_Thema4_ZitoumenoB

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