function highordersin

% Integration von sin im Intervall [0,pi]
% Achtung: Schlechte Implementierung!
a = 0;
b = pi;

err1 = [];
err7 = [];
err8 = [];
for n = 1:50,
    h = (b-a)/n;
    xi = a;
    q1 = 0; q7 = 0; q8 = 0;
    for i = 1:n,
       q1 = q1 + trapez(@sin, xi, xi+h);
       q7 = q7 + ncc7(@sin, xi, xi+h);
       q8 = q8 + ncc8(@sin, xi, xi+h);
       xi = xi + h;
    end
    err1 = [err1,abs(q1-2)];
    err7 = [err7,abs(q7-2)];
    err8 = [err8,abs(q8-2)];
end
subplot(1,2,1);
loglog(1:n,err1,'k-',1:n,err7,'b-',1:n,err8,'r-')
legend('Trapez', 'NC 7', 'NC 8');
makelarge
xlabel('#Teilintervalle');

subplot(1,2,2);

loglog(2:n+1,err1,'k-',(1:n)*7+1,err7,'b-',(1:n)*8+1,err8,'r-')
legend('Trapez', 'NC 7', 'NC 8');
makelarge
xlabel('#Funktionsauswertungen');


function q = trapez(f, a, b )

q = (b-a)*(feval(f,a)+feval(f,b))/2;

function q = ncc7(f, a, b )

alpha = [ 751/17280, 3577/17280, 49/640, 2989/17280, 2989/17280, ...
          49/640, 3577/17280, 751/17280 ];
     
x = linspace(a,b,8);
q = dot(alpha,feval(f,x)) * (b-a);

function q = ncc8(f, a, b )

alpha = [ 989/28350, 2944/14175, -464/14175, 5248/14175, -454/2835, ...
          5248/14175, -464/14175, 2944/14175, 989/28350 ];
     
x = linspace(a,b,9);
y = feval(f,x) * (b-a);
q = dot(alpha,feval(f,x)) * (b-a);