% Metoda kvadraticke interpolace 06

funkce='-2*sin(x) + sin(2*x) - 2*sin(3*x)/3';  % ucelova funkce        
a=1; b=3;                      % krajni hodnoty intervalu I=[a,b] 


r=0.5; s=1.5; t=3;    % pocatecni trojice bodu r,s,t
i=0;                  % promenna urcena k pocitani iteraci 
chyba=1;              % "mira presnosti" - jak presny bude vysledek 

fplot(funkce,[a,b],'r')                % graf funkce
line([2.356,2.356],[-3.185,-2.585])    % vyznaceni skutecneho minima
hold on;

disp(' '); % vlozi prazdny radek na vystupni obrazovku

% vypis na obrazovku

disp(sprintf('Je zadana funkce f=%s a body r=%f, s=%f a t=%f.',funkce,r,s,t));
disp(sprintf('Nasim cilem je najit minimum funkce f s presnosti 0,0001.'));

disp(' '); % vlozi prazdny radek na vystupni obrazovku

x=s;

% zakladni WHILE-cyklus (cyklus jednotlivych kroku) konci ve chvili
% dosazeni pozadovane presnosti, ci pokud citac iteraci dosahne hodnoty 15

while (abs(chyba>0.0001) & (i<15))
    
   i=i+1;  % zvyseni poctu iteraci     
   xx=x;   % ulozeni "stare" aproximace
   
   x=r;
   fr=eval(funkce);     % funkcni hodnota v bode r
   x=s;
   fs=eval(funkce);     % funkcni hodnota v bode s
   x=t;
   ft=eval(funkce);     % funkcni hodnota v bode t
  

   % vypocet minima prokladane paraboly 
   % promenne citatel a jmenovatel zavedeny kvuli prehlednosti
   
   citatel = fr*(t^2 - s^2) + fs*(r^2 - t^2) + ft*(s^2 - r^2);
   jmenovatel = fr*(t - s) + fs*(r - t) + ft*(s - r);
   
   % aby nedoslo k deleni nulou (nejen pri zaokrouhlovani by se mohlo stat)
   if (jmenovatel~=0) m=(1/2)*(citatel/jmenovatel);
   end    
   
   podminka=0;    % pomocna promenna, ktera hlida, jestli m==s
      
   x=m;
   fm=eval(funkce);     % funkcni hodnota v bode m
   
   disp(sprintf('%d. krok: f(r)=%f, f(s)=%f, f(t)=%f =>',i,fr,fs,ft));
   disp(sprintf('        => m=%f, f(m)=%f => ',m,fm));
   
   
   % vysetreni 3 moznosti, kam muze "padnout" bod m
     
                               
                                       
   if (m<s)
               if (fm<fs) % r=r;
                            t=s;         
                            s=m;                           
                           
               else         r=m; 
                          % s=s;        
                            t=t;
                           
               end    
              
       elseif (m==s) podminka=1;
           
       else %(m>s)
               if (fm<fs)   r=s; 
                            s=m;      
                          % t=t;
                           
               else       % r=r;
                          % s=s;     
                            t=m;
                           
               end            
               
                      
   end % konec podmineneho prikazu if  (radky 61-96) 
   
  
   % vyznaceni aproximace minima  
   
 line([x,x],[eval(funkce)-(0.02*(i+2)),eval(funkce)+(0.02*(i+2))])   

  % jestlize m==s, while-cyklus musi skoncit, "neni co dal pocitat"
      
  if (podminka)    i=17;                   
  end
        
  
  chyba = abs(xx - x);     % chyba jako rozdil poslednich dvou aproximaci  
  
  
  % vypis informaci o danem kroku na obrazovku
    
  disp(sprintf('        => r=%f, s=%f, t=%f, x=%f, chyba=%f',r,s,t,x,chyba)); 
  
  disp(' '); % vlozi prazdny radek na vystupni obrazovku
  keyboard;
end % konec while-cyklu (cyklu kroku vypoctu, radky 27-123)


% Vypis informaci o vypocitanem reseni v zavislosti na tom,
% jestli nastalo v prubehu vypoctu m==s

if (podminka)
    x=s-0.0001;
    fs1=eval(funkce);
    x=s+0.0001;
    fs2=eval(funkce);
    
    if ((fs<=fs1)&(fs<=fs2))     % spocitane minimum je skutecne min
                                 % az na danou presnost
       x=s;       
       disp(sprintf('Priblizne reseni: x=%f s presnosti 0.0001.',x));
       disp(sprintf('Hodnota minima: f(x)=%f.' ,eval(funkce)));

       
    else   % bod s byl "jen" minimem prokladane paraboly, ne uz funkce f
        
        disp(sprintf('Bod s byl nevhodne zvoleny. Je treba zvolit jiny.'));
    end
    
else  % m==s behem vypoctu nenastalo
    
disp(sprintf('Priblizne reseni: x=%f.',x));
disp(sprintf('Hodnota minima: f(x)=%f. ',eval(funkce)));
   
end