% Metoda regula falsi 06

funkce='-2*sin(x) + sin(2*x) - 2*sin(3*x)/3';    % ucelova funkce        
DER1='-2*cos(x) + 2*cos(2*x) - 2*cos(3*x)';      % derivace ucelove funkce

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

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

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

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

x=2.5;
fx1=eval(funkce)-0.08;
fx2=eval(funkce)+0.08;
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
   
   chyba = abs(x_k-x);     % chyba jako rozdil poslednich dvou aproximaci
   
   % 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, chyba=%f',x,chyba));
   end       
  
  
   disp(' '); % vlozi prazdny radek na vystupni obrazovku
  
   fx1=eval(funkce)-0.04*(i+2);
   fx2=eval(funkce)+0.04*(i+2);
   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