Algorithmensammlung: Numerik: Quadratur: Adaptive Multilevel-Quadratur
Adaptive Multilevel-Quadratur
BearbeitenDieses Quadraturverfahren basiert auf der summierten Trapezregel. Es verfeinert das Gitter entsprechend einer auf der Simpsonregel basierenden Fehlerabschätzung selbstständig.
Matlab
BearbeitenSiehe Octave.
Octave
Bearbeitenfunction [gitter, integral] = adaptquad(f, a, b, maxFehler, maxSchritte)
% Adaptive Multilevel-Quadratur basierend auf der Trapezregel
%
% f: Eine auf [a,b] definierte Funktion
% a,b: Grenzen
% maxFehler: Maximal zu akzeptierender Fehler
% maxSchritte: Maximale Schrittzahl (Anzahl Verfeinerungen)
%
% Die beiden folgenden Funktionen geben Arrays zurück, die die Quadratur
% der inneren Intervalle enthalten.
% Trapezregel
T = @(gitter, f) arrayfun(@(k) 1/2*(f(gitter(k-1)) + f(gitter(k)))*(gitter(k) - gitter(k-1)), 2:length(gitter));
% Simpsonregel
S = @(gitter, f) arrayfun(@(k) (gitter(k) - gitter(k-1))/6 * (f(gitter(k-1)) + 4*f(0.5*(gitter(k) + gitter(k-1))) + f(gitter(k))), 2:length(gitter));
% Fange mit einem simplen Gitter an
gitter = [a b];
schritt = 0;
while(true)
% Quadratur bestimmen
trap = T(gitter, f);
simp = S(gitter, f);
fehler = abs(simp - trap);
absFehler = abs(sum(simp) - sum(trap));
if(absFehler < maxFehler)
break
end;
if exist('maxSchritte') && schritt >= maxSchritte
break
end
% Gitter da wo der Fehler am größten ist verfeinern
newgitter = [];
for index = 1:length(gitter)
if (index > 1 && fehler(index-1) == max(fehler))
newgitter = [ newgitter (1/2*(gitter(index) + gitter(index-1))) ];
end;
newgitter = [ newgitter gitter(index) ];
end
gitter = newgitter;
schritt = schritt + 1;
end;
integral = sum(trap);
Quellen
BearbeitenDie Theorie zu diesem Algorithmus ist dem folgenden Vorlesungsscript entnommen:
Ralf Kornhuber, Christof Schütte; AG Numerische Mathematik (Hrsg.): Einführung in die Numerische Mathematik. Freie Universität Berlin April 2008, PDF, 2.4 MB