Zum Inhalt springen

Algorithmensammlung: Numerik: Hermiteinterpolation

Aus Wikibooks

Algorithmensammlung: Numerik

Hermiteinterpolation

[Bearbeiten]

Die Hermiteinterpolation ist ein Interpolationsverfahren, das in der Lage ist, Ableitungen der zu Interpolierenden Funktion zu berücksichtigen.

Der Algorithmus basiert auf der  klassischen Polynominterpolation und erweitert die dividierten Differenzen mithilfe der Hermite-Genocchi-Formel. In dem Fall, dass keine zwei Stützstellen zusammenfallen, ist das Ergebnis zu den Koeffizienten, die man von den dividieren Differenzen erhält, äquivalent. Für die mathematischen Hintergründe verweisen wir auf  Hermiteinterpolation.

Die Auswertung der ermittelten Stützstellen erfolgt mithilfe des Horner-Schemas.

Pseudocode

[Bearbeiten]
xvals ← Stützstellen
yvals ← Funktionswerte f(x) und ggf. Ableitungen bei mehrfachen x-Werten
zvals ← { f(xvals[i]) | i ∈ 1..#xvals }

for i ← #xvals..1 do
  for ji..#xvals do
    if i = j then
      [xi..xj]fzvals[i]
    else if xvals[i] = xvals[j] then
      index ← Index des ersten Vorkommens von xvals[i] in xvals
      [xi..xj]fyvals[j - i + index] / (j-i)!
    else
      [xi..xj]f ← ([xi+1..xj]f - [xi..xj-1]f) / (xvals[j] - xvals[i])
#include <stdlib.h>

void divdiff(double *x, double *y, double *stuetz, int n) {
    int i, j, k, l, m;
    double *koeffSchema;

    // koeffSchema repräsentiert eine n x n Matrix, zeilenweise
    // durchnummeriert.
    koeffSchema = malloc(sizeof(double) * (n * n));
    if(!koeffSchema) return;

    for(i=n-1; i>=0; i--) {
        for(j=i; j<n; j++) {
            if(i == j) {
                // Ersten Index finden, an dem x[i] vorkommt
                k = 0; while(x[i] != x[k]) k++;
                koeffSchema[i*n + j] = y[k];
            }
            else if(x[i] == x[j]) {
                // Fakultät von j-i bestimmen
                l = 1; m = j - i;
                while(m > 0) l *= m--;
                // Wiederum den ersten Index finden
                k = 0; while(x[i] != x[k]) k++;
                koeffSchema[i*n + j] = y[j - i + k] / l;
            }
            else {
                // Das ist wie im klassischen Fall
                koeffSchema[i*n + j] = (koeffSchema[i*n + n + j] - koeffSchema[i*n + j - 1]) /
                    (x[j] - x[i]);
            }
        }
    }

    for(i=0; i<n; i++) {
        stuetz[i] = koeffSchema[i];
    }

    free(koeffSchema);
}
public static double[] divdiff(double[] x, double[] y){
		double[] stuetz = new double [ y.length - 1 ];
		stuetz = y;
		int n = stuetz.length - 1;
		for(int i = n; i > 0; i--){
			for (int j = n; j > n - i ; j--){
				stuetz [ j ] = (stuetz [ j ] - stuetz [j - 1]) / (x [ j ] - x [ j-(n-i+1) ]);
			}
			
		}
		return(stuetz);
	}

Matlab

[Bearbeiten]

Siehe Octave.

Octave

[Bearbeiten]
function stuetzstellen = divdiff(x, y)
    n = length(x);
    % Aus dem Funktionswertevektor die jeweils ersten Werte zu einem bestimmten x finden
    z = arrayfun(@(i) y(min(find(x == x(i)))), 1:n);
    divdiff = zeros(n);
    for i = n:-1:1
        for j = i:n
            if i == j
                divdiff(i,j) = z(i);
            elseif x(i) == x(j)
                % Dies ist der Sonderfall in der Hermiteinterpolation
                anfangsIndex = min(find(x == x(i)));
                funktionsWert = y(j - i + anfangsIndex);
                faktor = factorial(j - i);
                divdiff(i, j) = funktionsWert / faktor;
            else
                divdiff(i, j) = (divdiff(i + 1, j) - divdiff(i, j - 1)) / (x(j) - x(i));
            end;
        end;
    end;
    stuetzstellen = divdiff(1,:);
end;