1D-Polynom-Interpolation

Satz (Polynom-Interpolation): Seien \((x_i, y_i) \in \real ^2\) für \(i = 0, \dotsc , n\) gegeben, wobei \(x_i \not = x_j\) für \(i \not = j\). Dann gibt es genau ein interpolierendes Polynom \(P_n(x) = \sum _{j=0}^n a_j x^j\) vom Grad \(\le n\).

Vandermonde-Matrix: Die Koeffizienten \(a_0, \dotsc , a_n\) lassen sich als Lösung des LGS
\(\smallpmatrix {1&x_0&\dots &x_0^n\\\vdots &\vdots &&\vdots \\1&x_n&\dots &x_n^n} \smallpmatrix {a_0\\\vdots \\a_n} = \smallpmatrix {y_0\\\vdots \\y_n} \iff V\vec {a} = \vec {y}\) bestimmen. \(V\) heißt Vandermonde-Matrix. Allerdings ist \(V\) schwierig und teuer zu berechnen und weitere Interpolationspkt.e sind aufwendig.

Lagrange-Interpolation: Die Lagrange-Interpolation bestimmt \(P_n\) durch den Ansatz \(P_n(x) = \sum _{i=0}^n y_i L_i(x)\) mit den Lagrange-Polynomen \(L_i(x) := \prod _{j\not =i} \frac {x - x_j}{x_i - x_j}\). \(L_i\) ist ein Polynom vom Grad \(n\) mit \(L_i(x_k) = \delta _{i,k}\).

Newton-Interpolation: Die Newton-Interpolation bestimmt \(P_n\) durch den Ansatz
\(P_n(x) = \sum _{i=0}^n a_i N_i(x)\) mit den Newton-Polynomen \(N_0(x) := 1\) und \(N_i(x) := \prod _{j=0}^{i-1} (x - x_j)\) für \(i = 1, \dotsc , n\). Wegen der rekursiven Definition (\(N_i(x) = (x - x_{i-1}) N_{i-1}(x)\)) gilt \(N_i(x_k) = 0\) für \(i > k\). Damit ist das resultierende LGS für die \(a_i\) in unterem Dreiecksformat und ein zusätzlicher Samplepunkt \(x_{n+1}\) verändert \(a_0, \dotsc , a_n\) nicht. Außerdem ist dieser Ansatz numerisch stabiler.

Kubische 1D-Interpolation

kubische 1D-Interpolation: Seien \(y_0, y_0’, y_1, y_1’ \in \real \) gegeben. Gesucht ist ein höchstens kubisches Polynom \(f(x) = ax^3 + bx^2 + cx + d\) mit \(f(i) = y_i\) und \(f’(i) = y_i’\) für \(i = 0, 1\). Man erhält das reguläre LGS \(A (a, b, c, d)^\tp = \vec {v}\) mit \(A := \smallpmatrix {0&0&0&1\\0&0&1&0\\3&2&1&0\\1&1&1&1}\) und \(\vec {v} := (y_0, y_0’, y_1’, y_1)^\tp \). Damit bekommt man \(f(x) = (a, b, c, d) (x^3, x^2, x, 1)^\tp = \vec {v}^\tp A^{-\tp } (x^3, x^2, x, 1)^\tp \) mit den
kubischen Hermite-Polynomen \(\smallpmatrix {H_0^3(x)\\H_1^3(x)\\H_2^3(x)\\H_3^3(x)} = A^{-\tp } (x^3, x^2, x, 1)^\tp = \smallpmatrix {(2x+1)(1-x)^2\\x(1-x)^2\\-x^2(1-x)\\(3-2x)x^2}\)

Bikubische Interpolation

bikubische Interpolation: Gegeben seien \(\vecs {f}{i,j} = (x_i, y_j)^\tp \) für \(x_i, y_j \in \{0, 1\}\) und Werte
\(z_{i,j}, \partial _x z_{i,j}, \partial _y z_{i,j}, \partial _x \partial _y z_{i,j} \in \real \). Gesucht ist bei der bikubischen Interpolation das bikubische Polynom \(f(x, y) = \sum _{n=0}^3 \sum _{m=0}^3 a_{n,m} x^n y^m\) mit \(f(\vecs {p}{i,j}) = z_{i,j}\), \(\partial _x f(\vecs {p}{i,j}) = \partial _x z_{i,j}\), \(\partial _y f(\vecs {p}{i,j}) = \partial _y z_{i,j}\) und \(\partial _x \partial _y f(\vecs {p}{i,j}) = \partial _x \partial _y z_{i,j}\). Schreibt man \(\vec {\alpha } := (a_{0,0}, a_{1,0}, a_{2,0}, a_{3,0}, a_{0,1}, \dotsc , a_{3,3})^\tp \) und
\(\vec {z} = (z_{0,0}, z_{1,0}, z_{0,1}, z_{1,1}, \partial _x z_{0,0}, \dotsc , \partial _x \partial _y z_{1,1})^\tp \), dann erhält man ein LGS \(B\vec {\alpha } = \vec {z}\) mit einer \((16 \times 16)\)-Matrix \(B\). Es gilt dann \(\vec {\alpha } = B^{-1} \vec {z}\).

bikubische Interpolation mit Hermite-Polynomen: Man kann den Ansatz auch schreiben als \(f(x, y) = \vec {y}^\tp A \vec {x}\) mit \(\vec {x} := (1, x, x^2, x^3)^\tp = C^\tp \vecs {h}{x}\), \(\vecs {h}{x} := (H_0^3(x), \dotsc , H_3^3(x))^\tp \) und einer bestimmten Basiswechsel-Matrix \(C = \smallpmatrix {1&0&0&0\\0&1&0&0\\0&1&2&3\\1&1&1&1}\), die man aus \(C^{-1} = \smallpmatrix {1&0&0&0\\0&1&0&0\\-3&-2&-1&3\\2&1&1&-2}\) erhält (analog \(\vec {y} = C^\tp \vecs {h}{y}\)). Damit bekommt man \(f(x, y) = \vecs {h}{y}^\tp C A C^\tp \vecs {h}{x}\). Man kann zeigen, dass
\(CAC^\tp = F := \smallpmatrix {f(0,0)&f_x(0,0)&f_x(1,0)&f(1,0)\\ f_y(0,0)&f_{xy}(0,0)&f_{xy}(1,0)&f_y(1,0)\\ f_y(0,1)&f_{xy}(0,1)&f_{xy}(1,1)&f_y(1,1)\\ f(0,1)&f_x(0,1)&f_x(1,1)&f(1,1)}\) bzw. \(A = C^{-1} F C^{-\tp }\).

Interpolation auf Dreiecken

baryzentrische Koordinaten: Seien \(\vec {a}, \vec {b}, \vec {c} \in \real ^2\) nicht kollinear. Dann gibt es für \(\vec {p} \in \real ^2\) baryzentrische Koordinaten \(\alpha _1, \alpha _2, \alpha _3 \in \real \) mit \(\vec {p} = \alpha _1 \vec {a} + \alpha _2 \vec {b} + \alpha _3 \vec {c}\) und \(\alpha _1 + \alpha _2 + \alpha _3 = 1\).
Durch Umschreiben \(\vec {p} - \vec {c} = \alpha _1 (\vec {a} - \vec {c}) + \alpha _2 (\vec {b} - \vec {c})\) erhält man \((\alpha _1, \alpha _2)^\tp = T^{-1} (\vec {p} - \vec {c})\) mit \(T := \smallpmatrix {\vec {a} - \vec {c} & \vec {b} - \vec {c}}\). Durch Ausschreiben der Inverse und Multiplikation bekommt man
\(\alpha _1 = \frac {1}{D} \det (\smallpmatrix {\vec {p} - \vec {c} & \vec {b} - \vec {c}})\), \(\alpha _2 = \frac {1}{D} \det (\smallpmatrix {\vec {a} - \vec {c} & \vec {p} - \vec {c}})\) und \(\alpha _3 = \frac {1}{D} \det (\smallpmatrix {\vec {b} - \vec {a} & \vec {p} - \vec {a}})\) mit der Determinate \(D := \det (T)\) (doppelter Flächeninhalt des Dreiecks).

geometrische Interpretation: \(\alpha _1\) ist der Anteil der Fläche des Dreiecks mit der zu \(\vec {a}\) gegenüber liegender Kante \(\vec {b} - \vec {c}\) als Kante und Ecke \(\vec {p}\) am gesamten Dreieck (falls \(\vec {p}\) im Dreieck liegt). Ob \(\vec {p}\) im Dreieck, auf einer Ecke/Kante oder außerhalb liegt, kann man leicht an den baryzentrischen Koordinaten ablesen.

Transformation zu Einheitsdreieck: Mit der Transformation \(x = a_x + (b_x - a_x) \xi + (c_x - a_x)\eta \) und \(y = a_y + (b_y - a_y) \xi + (c_y - a_y)\eta \) transformiert man ein Dreieck mit den Ecken \(\vec {a}, \vec {b}, \vec {c}\) auf das Einheitsdreieck. Im Einheitsdreieck gilt \(\alpha _1 = 1 - \xi - \eta \), \(\alpha _2 = \xi \) und \(\alpha _3 = \eta \).

baryzentrische Interpolation: Seien Werte \(f_1, f_2, f_3 \in \real \) an den Ecken des Dreiecks gegeben. Dann erhält man baryzentrische Interpolation durch \(f(\vec {p}) = \sum _{i=1}^3 \alpha _i(\vec {p}) f_i\).

quadratische Interpolation: Für quadr. Interpolation im Standarddreieck sei der Ansatz \(\vec {f}(\xi , \eta ) = \sum _{i=1}^6 f_i N_i(\xi , \eta )\) gegeben mit \(N_i(\xi _j, \eta _j) = \delta _{i,j}\) für \((\xi _1, \eta _1) := (0, 0)\), \((\xi _2, \eta _2) := (1, 0)\), \((\xi _3, \eta _3) := (0, 1)\), \((\xi _4, \eta _4) := (1/2, 0)\), \((\xi _5, \eta _5) := (1/2, 1/2)\) und \((\xi _6, \eta _6) := (0, 1/2)\).
Man erhält \(N_1(\xi , \eta ) = (1 - \xi - \eta ) (1 - 2\xi - 2\eta )\), \(N_2(\xi , \eta ) = \xi (2\xi - 1)\), \(N_3(\xi , \eta ) = \eta (2\eta - 1)\), \(N_4(\xi , \eta ) = 4\xi (1 - \xi - \eta )\), \(N_5(\xi , \eta ) = 4\xi \eta \) und \(N_6(\xi , \eta ) = 4\eta (1 - \xi - \eta )\).

kubische Interpolation: Man kann den Ansatz auch für höhere Ordnungen verallgemeinern. Zum Beispiel wäre der kubische Ansatz
\(u(\xi , \eta ) = c_1 + c_2\xi + c_3\eta + c_4\xi ^2 + c_5\xi \eta + c_6\eta ^2 + c_7\xi ^3 + c_8\xi ^2\eta + c_9\xi \eta ^2 + c_{10}\eta ^3\).

Bikubische Interpolation auf krummlinigen Gittern

krummlinige Gitter: Zur Interpolation auf krummlinigen Gittern muss man die Ableitungen, die im physischen Raum (P-Raum) mit Koord.en \((x, y)\) gegeben sind, umrechnen auf den Berechnungsraum (C-Raum) mit Koord.en \((\xi , \eta )\). Im C-Raum sollte das Gitter ein uniformes 2D-Rechtecksgitter mit Gitterabständen \(\Delta \xi = 1 = \Delta \eta \) sein.

Umrechnung der Ableitungen: Für \(\vecs {\partial }{xy} := (\partial _x, \partial _y, \partial _x^2, \partial _x \partial _y, \partial _y^2)^\tp \) und
\(\vecs {\partial }{\xi \eta } := (\partial _\xi , \partial _\eta , \partial _\xi ^2, \partial _\xi \partial _\eta , \partial _\eta ^2)^\tp \) erhält man \(\vecs {\partial }{xy} = B\vecs {\partial }{\xi \eta }\) mit \(B := \smallpmatrix {B_1&0\\B_2&B_3}\) und
\(B_1 := \smallpmatrix {\partial _x \xi &\partial _x \eta \\\partial _y \xi &\partial _y \eta }\), \(B_2 := \smallpmatrix {\partial _x^2 \xi &\partial _x^2 \eta \\\partial _x \partial _y \xi & \partial _x \partial _y \eta \\\partial _y^2 \xi &\partial _y^2 \eta }\) sowie \(B_3 := \smallpmatrix {(\partial _x \xi )^2&2\partial _x \xi \cdot \partial _x \eta & (\partial _x \eta )^2\\ \partial _x \xi \cdot \partial _y \xi & \partial _x \xi \cdot \partial _y \eta + \partial _x \eta \cdot \partial _y \xi & \partial _x \eta \cdot \partial _y \eta \\ (\partial _y \xi )^2&2\partial _y \xi \cdot \partial _y \eta &(\partial _y \eta )^2}\).
Umgekehrt gilt \(\vecs {\partial }{\xi \eta } = C\vecs {\partial }{xy}\) mit \(C := B^{-1} = \smallpmatrix {C_1&0\\C_2&C_3}\) und
\(C_1 := B_1^{-1} = \smallpmatrix {\partial _\xi x&\partial _\xi y\\\partial _\eta x&\partial _\eta y}\), \(C_2 := -B_3^{-1} B_2 B_1^{-1} = \smallpmatrix {\partial _x^2 x&\partial _\xi ^2 y\\\partial _\xi \partial _\eta x& \partial _\xi \partial _\eta y\\\partial _\eta ^2 x&\partial _\eta ^2 y}\) sowie
\(C_3 := B_3^{-1} = \smallpmatrix {(\partial _\xi x)^2&2\partial _\xi x \cdot \partial _\xi y& (\partial _\xi y)^2\\ \partial _\xi x \cdot \partial _\eta x& \partial _\xi x \cdot \partial _\eta y + \partial _\xi y \cdot \partial _\eta x& \partial _\xi y \cdot \partial _\eta y\\ (\partial _\eta x)^2&2\partial _\eta x \cdot \partial _\eta y&(\partial _\eta y)^2}\).