function [vals,errs] = aufg03f01_08

%AUFG02F06 implementiert und testet verschiedene Quadraturformeln.
%
%   [VALS, ERRS] = AUFG03F01 implementiert am Beispiele der
%   drei Integrale
%
%           1
%           /
%      I1 = | (1-1/exp(x)) dx
%           /
%           0
%
%           1
%           /
%      I2 = | 2x/(1+x^2) dx
%           /
%           0
%
%           1
%           /
%      I3 = | sqrt(x) dx
%           /
%           0
%
%   folgende Quadraturformeln:
%
%     a) abgeschlossene Newton-Cotes Formeln fuer n=1..4
%     b) die summierte Simpson-Regel fuer m=2k+1, k=5..10
%     c) Legendre-Gauss-Quadratur fuer tabellierte Gewichte
%        und Knoten fuer n=1..5
%
%   VALS und ERRS sind Arrays der Dimensionen 15 x 3 und
%   enthalten die berechneten Werte und den Fehler dieser
%   Naeherungen, respektive. In einer Zeile stehen die Werte
%   fuer ein Verfahren, in den Spalten stehen die Naeherungen
%   fuer ein Integral. Sortiert sind die Felder nach der
%   Abarbeitungsreihenfolge (I1 bis I3, Newton-Cotes, summierte
%   Simpson-Regel, Legendre-Gauss, kleine Ordnung zuerst).
%
%   AUFG03F01_08 ist Teil der Loesung von Aufgabe 3, Blatt 3
%   der Vorlesung 'Numerische Verfahren' von Jens Zemke. Die
%   Funktion arbeitet unter Matlab und Octave.
%
%   SEE ALSO: quad, quadl

%   $Author: zemke $
%   $Revision: 1.3 $
%   $Date: 2003/04/22 13:15:16 $

%% for better readability ;-)
format long

vals = zeros(15,3);
errs = vals;

exact = [exp(-1) log(2) 2/3]


%% closed Newton-Cotes
% Insel-Mathematik :: Verwendbar und durchfuehrbar im Kopf
%

cnc = zeros(4,3);

%
% weights + nodes
%

w1 = 0.5 *[1 1];
w2 = 1/6 *[1 4 1];
w3 = 1/8 *[1 3 3 1];
w4 = 1/90*[7 32 12 32 7];

n1 = linspace(0,1,2)';
n2 = linspace(0,1,3)';
n3 = linspace(0,1,4)';
n4 = linspace(0,1,5)';

%
% closed Newton-Cotes approximations
%

cnc(1,1) = w1*(1-1./exp(n1));
cnc(2,1) = w2*(1-1./exp(n2));
cnc(3,1) = w3*(1-1./exp(n3));
cnc(4,1) = w4*(1-1./exp(n4));

cnc(1,2) = w1*((2*n1)./(1+n1.^2));
cnc(2,2) = w2*((2*n2)./(1+n2.^2));
cnc(3,2) = w3*((2*n3)./(1+n3.^2));
cnc(4,2) = w4*((2*n4)./(1+n4.^2));

cnc(1,3) = w1*sqrt(n1);
cnc(2,3) = w2*sqrt(n2);
cnc(3,3) = w3*sqrt(n3);
cnc(4,3) = w4*sqrt(n4);

vals(1:4,:) = cnc;
disp('Closed Newton-Cotes-Quadrature, n=1...4')
cnc

%% summed Simpson formula
% Rapid-Prototyping :: Schnell implementiert, zum Testen
%

i=1;
ssf = zeros(6,3);

for k = 5:10

%
% generate nodes => simple.
%
   nodes = linspace(0,1,2*k+1)'; %'
%
% generate weights => more or less simple.
%
   wss = zeros(1,2*k+1);
   wss([1 2*k+1])     = 1;
   wss(2*(1:k))       = 4;
   wss(2*(1:(k-1))+1) = 2;
   wss = 1/(6*k)*wss;
%
% evaluate => always as scalar-product.
%
   ssf(i,1) = wss*(1-1./exp(nodes));
   ssf(i,2) = wss*((2*nodes)./(1+nodes.^2));
   ssf(i,3) = wss*sqrt(nodes);
   i = i+1;

end;

vals(5:10,:) = ssf;
disp('Summed Simpson-Formula, m=2*k, k=5...10');
ssf


%% Gauss-quadrature
% Novell Cuisine ;-) :: Komplex, aber gut
%

%
% table of weights for [-1,1]
%
wg1 = 2;
wg2 = [1 1];
wg3 = 1/9*[5 8 5];
wg4 = 1/36*(18+sqrt(30)*[-1 1 1 -1]);
wg5 = 1/900*([322-13*sqrt(70) 322+13*sqrt(70) 512 ...
              322+13*sqrt(70) 322-13*sqrt(70) ]);

%
% table of nodes for [-1,1]
%
ng1 = 0;
ng2 = 1/sqrt(3)*[-1 1]';
ng3 = sqrt(3/5)*[-1 0 1]';
ng4 = 1/35*([-sqrt(525+70*sqrt(30)) -sqrt(525-70*sqrt(30)) ...
              sqrt(525-70*sqrt(30))  sqrt(525+70*sqrt(30))]');
ng5 = 1/21*([-sqrt(245+14*sqrt(70)) -sqrt(245-14*sqrt(70)) ...
            0 sqrt(245-14*sqrt(70))  sqrt(245+14*sqrt(70))]');

%
% transform nodes to [0,1]
%
ngt1 = 0.5*ng1+0.5;
ngt2 = 0.5*ng2+0.5;
ngt3 = 0.5*ng3+0.5;
ngt4 = 0.5*ng4+0.5;
ngt5 = 0.5*ng5+0.5;

%
% evaluate, transformation of weights included (0.5*wi)
%
lgq = zeros(5,3);

lgq(1,1) = 0.5*wg1*(1-1./exp(ngt1));
lgq(2,1) = 0.5*wg2*(1-1./exp(ngt2));
lgq(3,1) = 0.5*wg3*(1-1./exp(ngt3));
lgq(4,1) = 0.5*wg4*(1-1./exp(ngt4));
lgq(5,1) = 0.5*wg5*(1-1./exp(ngt5));

lgq(1,2) = 0.5*wg1*((2*ngt1)./(1+ngt1.^2));
lgq(2,2) = 0.5*wg2*((2*ngt2)./(1+ngt2.^2));
lgq(3,2) = 0.5*wg3*((2*ngt3)./(1+ngt3.^2));
lgq(4,2) = 0.5*wg4*((2*ngt4)./(1+ngt4.^2));
lgq(5,2) = 0.5*wg5*((2*ngt5)./(1+ngt5.^2));

lgq(1,3) = 0.5*wg1*sqrt(ngt1);
lgq(2,3) = 0.5*wg2*sqrt(ngt2);
lgq(3,3) = 0.5*wg3*sqrt(ngt3);
lgq(4,3) = 0.5*wg4*sqrt(ngt4);
lgq(5,3) = 0.5*wg5*sqrt(ngt5);

vals(11:15,:) = lgq;
disp('Legendre-Gauss-Quadrature, Table-driven, n=1...5');
lgq

%
% compute errors.
%
errs = vals-ones(15,1)*exact;
