Algorithmensammlung: Numerik: Hermiteinterpolation

Hermiteinterpolation

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

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 j ← i..#xvals do
    if i = j then
      [xi..xj]f ← zvals[i]
    else if xvals[i] = xvals[j] then
      index ← Index des ersten Vorkommens von xvals[i] in xvals
      [xi..xj]f ← yvals[j - i + index] / (j-i)!
    else
      [xi..xj]f ← ([xi+1..xj]f - [xi..xj-1]f) / (xvals[j] - xvals[i])

C

#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);
}

Java

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

Siehe Octave.

Octave

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;