Algorithmensammlung: Numerik: Hermiteinterpolation
Aus Wikibooks
-
- Dividierte Differenzen
- Hermiteinterpolation
- Horner-Schema
- Quadratur
[Bearbeiten] 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.
Inhaltsverzeichnis |
[Bearbeiten] 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 klassischem 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); }
[Bearbeiten] Matlab
Siehe Octave.
[Bearbeiten] 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;
[Bearbeiten] 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])