% Metoda regula falsi 04

funkce='x + 3/(x^2)';     % ucelova funkce        
DER1='1-6/(x^3)';         % derivace ucelove funkce

a=0.5; b=3;               % krajni hodnoty intervalu I=[a,b] 
x_s=1.75;                 % 1.aproximace 
x_k=2.375;                % 2.aproximace

% body x_k a x_s hraji roly stejnojmennych bodu ve vzorci
% pro vypocet (k+1)-ni aproximace

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.326,3.126])    % vyznaceni skutecneho minima

x=1.75;
fx1=eval(funkce)-0.1;
fx2=eval(funkce)+0.1;
line([x,x],[fx1,fx2])    % vyznaceni aproximace minima

x=2.375;
fx1=eval(funkce)-0.15;
fx2=eval(funkce)+0.15;
line([x,x],[fx1,fx2])    % vyznaceni aproximace minima


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

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


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

while ((chyba>0.0001) & (i<15))
    
   i=i+1;  % zvyseni poctu iteraci     
     
   disp(sprintf('%d. krok: x_k=%f, x_s=%f, ',i,x_k,x_s));
  
   x=x_s;
   DERx_s=eval(DER1);      % hodnota derivace funkce f v bode x_s
   x=x_k;   
   DERx_k=eval(DER1);      % hodnota derivace funkce f v bode x_k
  
   disp(sprintf('          df(x_k)=%f, df(x_s)=%f =>',DERx_k,DERx_s));
   
   % vzorec pro vypocet nasledujici iterace
   x = x_k - ((x_k-x_s)/(DERx_k-DERx_s))*(DERx_k);
   
   DERx=eval(DER1);      % hodnota derivace funkce f v bode x
   
   
   % bod x_s "opravime" (pokud je potreba) tak, 
   % aby v pristi iteraci platilo: f'(x_k)*f'(x_s)<0
   
   if (DERx*DERx_s>0) x_s=x_k;
   end
   
   % nenasli jsme uz nahodou minimum? - at ten vypis nevypada divne
   
   if (DERx_k==0)
     disp(sprintf('        => aproximace: x=%f je minimem',x));
   else
     disp(sprintf('        => nova aproximace: x=%f',x));
   end  
  
   chyba = abs(x_k-x);     % chyba jako rozdil poslednich dvou aproximaci
  
  
   disp(' '); % vlozi prazdny radek na vystupni obrazovku
  
   fx1=eval(funkce)-0.1*(i+1);
   fx2=eval(funkce)+0.1*(i+1);
   line([x,x],[fx1,fx2])    % vyznaceni aproximace minima
 
   
   if (chyba==0) i=17;
   end
  
   x_k=x; % ulozeni nove spocitane aproximace do x_k pro novy vypocet
   keyboard;
end % konec while-cyklu (cyklu kroku vypoctu, radky 29-74)


% Vypis informaci o vypocitanem reseni v zavislosti na tom,
% jestli jsme nahodou "netrefili" minimum presne

if (i==17)
      disp(sprintf('Presne reseni: x=%f.',x));
      disp(sprintf('Hodnota minima: f(x)=%f.' ,eval(funkce)));
      
else  
   
disp(sprintf('Priblizne reseni: x=%f s chybou %f.',x,chyba));
disp(sprintf('Hodnota minima: f(x)=%f. ',eval(funkce)));
   
end