% Metoda zlateho rezu 04

funkce='x + 3/(x^2)';                % ucelova funkce
a=0.5; b=3;                          % krajni hodnoty intervalu I=[a,b] 

p = linspace(a,b);                   % vykreslovany interval
q = (p.^(1)) + 3*(p.^(-2));          % predpis funkce
plot(p,q,'r')                        % vykresleni funkce   
line([1.817,1.817],[2.526,2.926])    % vyznaceni skutecneho minima

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

% vypis na obrazovku

disp(sprintf('Je zadana funkce f=%s a interval [a,b]=[%f,%f].',funkce,a,b));
disp(sprintf('Nasim cilem je najit minimum funkce f s presnosti 0,0001.'));

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

r=((5^(0.5))-1)/2;          % rez r

i=0;         % promenna urcena k pocitani iteraci
chyba=1;     % "mira presnosti" - jak presny bude vysledek


% cyklus skonci ve chvili, kdy bude dosazeno pozadovane presnosti
% nebo kdyz pocet kroku dosahne hodnoty 20 (oboji je mozne zvolit libovolne)

 while (chyba>0.0001) & (i<20)

    c = a + (1-r)*(b-a);       % body urcene rezem intervalu
    d = a + r*(b-a);           

    x=c; 
    fc = eval(funkce);         % funkcni hodnota v bode c
    x=d;
    fd = eval(funkce);         % funkcni hodnota v bode d
    
    if (fd<fc)  a=c;           % novy interval bude [c,b]
         else   b=d;           % novy interval bude [a,d]
    end                          


    % vyznaceni "smrstovani" puvodniho intervalu [a,b]
    % omezeni na i<20 je kvuli prehlednosti obrazku
    
    if (i<20)
        p = linspace(a,b);
        q2 = (p.^(1)) + 3*(p.^(-2)) - 0.3*(i+1);    
        hold on;
    end

    % podmineny prikaz zakresluje "smrstovani" dvema barvami (modra, cerna),
    % aby byly lepe rozlisene jednotlive kroky

    if (mod(i,2) == 0) plot(p,q2,'b')
            else       plot(p,q2,'k')
    end    


    chyba=(b-a)/2;            % chybu urcuje rozdil funkcnich hodnot
    i=i+1;                    % konci cyklus, zvysime citac iteraci


    % vypis informaci o krocich vypoctu na obrazovku

    disp(sprintf('%d. krok: c=%f, d=%f, f(c)=%f, f(d)=%f.',i,c,d,fc,fd));   
    disp(sprintf('          [a,b]=[%f, %f]. ',a,b));
    disp(' '); % vlozi prazdny radek na vystupni obrazovku
    keyboard;


end  % zde konci while-cyklus (cyklus kroku vypoctu)


x=(a+b)/2;  % konecna aproximace jako stred konecneho intervalu

% vypis konecneho vysledku na obrazovku

disp(sprintf('Priblizne reseni: x=%f s presnosti %f.',x,chyba));   
disp(sprintf('Hodnota minima: f(x)=%f. ',eval(funkce))); 
 
 
