% Integration von sqrt im Intervall [0,1]
a = 0;
b = 1;

nmax = 50;
err1 = zeros(1,nmax);
err2 = zeros(1,nmax);
err1g = zeros(1,nmax);
err2g = zeros(1,nmax);


for n = 1:nmax,
    h = (b-a)/n;
    xi = a;
    % Trapez
    q1 = h * ( (sqrt(a)+sqrt(b))/2 + sum( sqrt(a+(1:n-1)*h) ) );
    
    % Simpson
    f = sqrt(linspace(a,b,2*n+1));
    f(2:2:end-1) = 4*f(2:2:end-1);
    f(3:2:end-2) = 2*f(3:2:end-2);
    q2 = h * sum(f) / 6;

    % Gauss, n = 1
    hp = (1/2 - 1/(2*sqrt(3)))*h;
    x = a+(0:n-1)*h;
    q1g = h*( sum(sqrt(x+hp)) + sum(sqrt(x+h-hp)) ) / 2;

    % Gauss, n = 2
    hp = (1/2 - 1/2*sqrt(3/5))*h;
    x = a+(0:n-1)*h;
    q2g = 5*sum(sqrt(x+hp)) + 5*sum(sqrt(x+h-hp)) + 8*sum(sqrt(x+h/2));
    q2g = q2g * h / 18;


    err1(n) = abs(q1-2/3);
    err2(n) = abs(q2-2/3);
    err1g(n) = abs(q1g-2/3);
    err2g(n) = abs(q2g-2/3);
end
subplot(1,2,1);
loglog(1:n,err1,'k-',1:n,err2,'b-',1:n,err1g,'r-',1:n,err2g,'g-')
legend('Trapez', 'Simpson', 'Gauss, n = 1', 'Gauss, n = 2');
makelarge
xlabel('#Teilintervalle');

subplot(1,2,2);

loglog(2:n+1,err1,'k-',(1:n)*2+1,err2,'b-',(1:n)*2,err1g,'r-',(1:n)*3,err2g,'g-')
legend('Trapez', 'Simpson', 'Gauss, n = 1', 'Gauss, n = 2');
makelarge
xlabel('#Funktionsauswertungen');