%% Annnnna

close all % Close all open figure windows.

xi = [-3/2 -1/2 1/2 3/2]; 
f=@(x)x.^xi; % The function we evaluate.

a=1; b=2; 
n=10;
h=(b-a)/n;
dx = h; % alias


%% Left Riemann sum
% x=linspace(a,b,n+1);
% qv=sum(f(x(1:n)).*dx)

% subplot(2,2,1)
% plot(x, f(x), 'r-'), hold on, grid on
% for i=1:n
%     plot([x(i) x(i+1)], [f(x(i)) f(x(i))], 'g-')
%     plot([x(i+1) x(i+1)], [f(x(i)) f(x(i+1))], 'g-')
% end
% title('v�nster rektangelregel');
%     

%% Right Riemann sum
% x=linspace(a,b,n+1);
% qh=sum(f(x(2:n+1)).*dx)
% 
% subplot(2,2,2)
% plot(x, f(x), 'r-'), hold on, grid on
% for i=1:n
%     plot([x(i) x(i+1)], [f(x(i+1)) f(x(i+1))], 'g-')
%     plot([x(i) x(i)], [f(x(i)) f(x(i+1))], 'g-')
% end
% title('H�ger rektangelregel');

a=0; b=2;
xi = [-3/2 -1/2 1/2 3/2]; 
for blep = xi
    i = 1;
    for n=[10 20 40 80]
        x=linspace(a,b,n+1);
        h=(b-a)/n;
        f = @(x) x.^blep;
        kontroll = integral(f, a, b); % Accurate control integral for error computation.
        qm=sum(f((x(1:n)+x(2:n+1))./2).*h);
        err(i) = abs(qm - kontroll)
       
        i = i+1;
    end
    plot([10 20 40 80], err)
    hold on
    grid on
end
xlabel('n');ylabel('err(h)');

%% Middle Riemann sum
x=linspace(a,b,n+1);
qm=sum(f((x(1:n)+x(2:n+1))./2).*dx)
err = abs(qm - kontroll)

%subplot(2,2,3:4)
plot(x, f(x), 'r-'), hold on, grid on
for i=1:n
    plot([x(i) x(i)], [f(x(i)) (f(x(i))+f(x(i+1)))/2], 'g-')
    plot([x(i) x(i+1)], [(f(x(i))+f(x(i+1)))/2  (f(x(i))+f(x(i+1)))/2], 'g-')
    plot([x(i+1) x(i+1)], [(f(x(i))+f(x(i+1)))/2 f(x(i+1))], 'g-')
end

title(['Mittpunktsregeln rektangelregel. n= ', num2str(n), ' err= ', num2str(err)]);

%% Trapezoidal sum
x=linspace(a,b,n+1);
qt=sum(((f(x(1:n))+f(x(2:n+1)))./2).*dx);
%title('Trapetzmetoden')