Algorithmensammlung: Numerik: Hermiteinterpolation

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

Siehe 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;