% Integration von sin im Intervall [0,pi]
a = 0;
b = pi;

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 * ( (sin(a)+sin(b))/2 + sum( sin(a+(1:n-1)*h) ) );
    
    % Simpson
    f = sin(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(sin(x+hp)) + sum(sin(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(sin(x+hp)) + 5*sum(sin(x+h-hp)) + 8*sum(sin(x+h/2));
    q2g = q2g * h / 18;


    err1(n) = abs(q1-2);
    err2(n) = abs(q2-2);
    err1g(n) = abs(q1g-2);
    err2g(n) = abs(q2g-2);
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');