% CVICENI - 5

% Numericka metoda puleni intervalu
% funkce jsou ulozeny ve stejnem adresari jako CV5.m
x = Puleni("f1", -10,10)

% Numericka derivace:
% --------------------------------------%
% df/dx(x) = lim[ ( f(x+dx) - f(x) ) / dx ], dx-> 0
% Numericky: dx -> tolerance (nemuze byt 0)

y = df1("f",1,1e-07)

% ------------------------------------------------%
Numericka integrace:
Iexact = -5/6
a = 0;
b = 1;
N = 50;
x = linspace(a,b,N);
I = 0;

for i=1:N-1
    dx = x(i+1) - x(i);
    dF =  f(x(i)) * dx; 
    I  = I + dF;
end
I0 = I;

I = 0;
for i=2:N-1
    dx = x(i+1) - x(i);
    dF =  f(x(i)) * dx; 
    I  = I + dF;
end
I = I + (f(x(1)) + f(x(N)))/2 * dx;

I0
I1 = I
Iexact


% Graphics:
% ------------------------------------------------%
x_fine = linspace(a,b,1000);
eps = 1e-5;
Fig = figure(1);

x_fine(end)

for N=1:10:100    
    x = linspace(a,b,N);
    I = 0; 
    for i=1:N-1
        hold on;
        dx = x(i+1) - x(i);
        dF =  f(x(i)) * dx; 
        I  = I + dF;
        x0 = x(i);
        Dx = x(i+1)-x(i);
        if f(x(i)) >= 0
            dy = f(x(i));
            y0 = 0;
        else
            y0 = f(x(i));
            dy = abs(f(x(i)));
        end
        r = rectangle('Position', [x0,y0,dx,dy+eps], 'FaceColor', [0, 0, 1, 0.4], 'EdgeColor', 'k'); 
    end
    plot(x_fine, f(x_fine), "-r", "LineWidth", 4);
    value = sprintf('I =  %.5f',I);
    legend(value)
    hold off;
    pause(2);
    clf(Fig);
end

%%%-------------------------FUNKCE-------------------------%%%
function y = f(x)
    y = x.^5 - 3*x.^2;
end

function y = df1(f,x,dx)
    y =  (feval("f",x+dx) - feval("f",x)) / dx;
end
%---------------------------------------------------%
function y = f1(x)
    y = 10*x^5 - 3*x^3 + 2*x^2 + x + 5;
end
%---------------------------------------------------%
function [x]=Puleni(jmenof,a,b)
% Obecne reseni rovnice f(x)=0 metodou puleni intervalu
% [x]=Puleni(jmenof,a,b);
% x ... nalezene reseni rovnice
% jmenof ... nazev funkce f(x) jako retezec znaku
% a ... dolni mez intervalu
% b ... horni mez intervalu
nmax=40; 	% maximalni pocet vycisleni
fa=feval(jmenof,a); 	% levy okraj
fb=feval(jmenof,b); 	% pravy okraj
if fa*fb>0
   x=inf; 	% chybny predpoklad uspechu
   return;
end
for k=3:nmax 	% bezpecny cykl
   x=(a+b)/2; 	% stred intervalu
   f=feval(jmenof,x); 	% vycislit v polovine
   if fa*f<=0 	% je neco v intervalu <a;x>?
      b=x; 	% zahozeni intervalu (x;b>
   else
      a=x; 	% zahozeni intervalu <a;x)
   end
end
x=(a+b)/2; 	% naposled
end
