Algorithmensammlung: Numerik: Quadratur: Adaptive Multilevel-Quadratur

Algorithmensammlung: Numerik: Quadratur

• Trapezregel
• Simpson-Regel
• Romberg-Verfahren
• Newton-Cotes-Quadratur
• Adaptive Multilevel-Quadratur

Adaptive Multilevel-Quadratur

Dieses Quadraturverfahren basiert auf der summierten  Trapezregel. Es verfeinert das Gitter entsprechend einer auf der  Simpsonregel basierenden Fehlerabschätzung selbstständig.

Matlab

Siehe Octave.

Octave

function [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

Die 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