Numerik Interpolation PDF

Kapitel 2 Interpolation 2.1 Polynomiale Interpolation Im den ersten beiden Kapiteln dieser Vorlesung behandeln wir die Probleme der numerischen Inter- polation und Integration im univariaten bzw. eindimensionalen Fall (Betre↵s der Erweiterung auf den multivariaten bzw. mehrdimensionalen Fall w...

Kapitel 2 Interpolation 2.1 Polynomiale Interpolation Im den ersten beiden Kapiteln dieser Vorlesung behandeln wir die Probleme der numerischen Inter- polation und Integration im univariaten bzw. eindimensionalen Fall (Betre↵s der Erweiterung auf den multivariaten bzw. mehrdimensionalen Fall wird – bis auf vereinzelte Ausnahmen in den Übungen – auf weiterführende Veranstaltungen verwiesen). Eindimensionale Interpolationsaufgaben behandeln folgende Fragestellung: Gegeben sind n + 1 Paare (xi , yi ) 2 R ⇥ R, i = 0, 1,..., n mit xi 6= xk , i 6= k und eine Klasse reellwertiger Funktionen (·; c0 ,..., cn ), die von n + 1 Parametern c0 ,..., cn abhängen. Gesucht sind die Parameterwerte c0 ,..., cn 2 R mit der (Interpolations-) Eigenschaft (xi ; c0 ,..., cn ) = yi , i = 0, 1,..., n. (2.1) Man kann sich unter den Punktpaaren die Ergebnisse einer Messreihe an den Messstellen xi mit den Messwerten yi vorstellen. Das Problem (2.1) entspricht der Lösung eines i.a. Falle nichtlinearen Gleichungssystems. Speziell erhält man ein lineares System, wenn die Parameter c0 ,..., cn nur linear in eingehen. Mit einer geeigneten Basis { 0 ,... , n } gelte dann n X (x; c0 ,..., cn ) := ci i (x). (2.2) i=0 Wir gehen hier nur auf diesen linearen Fall ein. Aus historischer Sicht ist die polynomiale Interpolation der wichtigste Spezialfall. In der Monom-Basis mit i (x) = xi ergibt sich zum Beispiel n X (x; c0 ,..., cn ) := c i xi. i=0 Zunächst wurde diese Interpolationsform vorwiegend benutzt, um in Tafelwerken angegebene diskrete Funktionswerte zu interpolieren, vor allem durch lineare und quadratische Funktionen. Im Zeitalter elektronischer Rechentechnik hat das nur noch geringe Bedeutung, da die Funktionswerte durch schnelle Routinen verfügbar sind. Hingegen ist die polynomiale Interpolation nach wie vor wesentlich zur Begründung anderer numerischer Verfahren (z.B. numerische Integration, vgl. Kapitel 3.1 und 13 14 KAPITEL 2. INTERPOLATION 3.2, Di↵erentiation und Näherungslösung von Di↵erentialgleichungen oder allgemeiner von Opera- torgleichungen). Bei der Interpolation von Messreihen (xi , yi ), i = 0,..., n ist die polynomiale Interpolation vor allem bei größeren Werten von n wenig geeignet, da das interpolierende Polynom oft an den Intervallen- den stark oszilliert (vgl. Beispiel 2.4). Abhilfe scha↵t hier die Interpolation durch Splines, d.h. die Annäherung durch stückweise polynomiale Funktionen (vgl. Kapitel 2.3). Trigonometrische Polyno- me spielen bei der Interpolation periodischer Daten eine zentrale Rolle (vgl. Kapitel 2.2). Bei den einzelnen Interpolationsansätzen (vgl. Kapitel 2.1–??) behandeln wir Existenz- und Ein- deutigkeitsaussagen, dann die praktische Darstellung und Berechnung sowie teilweise Fehler- und Konvergenzaussagen. 2.1.1 Lagrangesche Interpolation Sei ⇧n der lineare Raum der Polynome vom Grad kleiner oder gleich n 2 N0 als Unterraum des linearen Raumes C[a, b] der stetigen Funktionen auf [a, b] mit a < b. Wir werden mehrfach Gebrauch von folgender Eigenschaft algebraischer Polynome machen. Satz 2.1 (i) Jedes Polynom aus ⇧n , n 2 N0 , das mehr als n (komplexe) Nullstellen (bei Berücksichtigung ihrer Vielfachheit) besitzt, verschwindet identisch. (ii) Die Monome i (x) = xi , i = 0,..., n bilden eine Basis des Raumes ⇧n , d.h. dim(⇧n ) = n + 1, ⇧n = span{1, x,..., xn }. Beweis. Aussage (i) beweist man, ausgehend vom Fundamentalsatz der Algebra, am einfachsten per Induktion über n. Aussage (ii) ist dann Folgerung von Aussage (i). (Übungsaufgabe !) Definition 2.2 Die Lagrangesche Interpolationsaufgabe lautet: Bestimme ein Polynom Pn 2 ⇧n zu gegebenen Stützstellen x0 ,..., xn 2 R mit xi 6= xk , i 6= k und Funktionswerten y0 ,..., yn , so dass Pn (xi ) = yi , i = 0,..., n. (2.3) Satz 2.3 Die Lagrangesche Interpolationsaufgabe ist eindeutig lösbar. Pn hat die Lagrange-Darstellung n X Pn (x) = yk lk (x) (2.4) k=0 mit den Lagrange-Polynomen n Y x xj lk (x) := , k = 0,..., n. (2.5) j =0 xk xj j 6= k 2.1. POLYNOMIALE INTERPOLATION 15 Beweis. Existenz: Per Ansatz gilt Pn 2 ⇧n sowie ⇢ 0, k 6= i lk (xi ) = ik =. 1, k = i Damit ist die Interpolationsbedingung erfüllt. Eindeutigkeit: Seien p1 , p2 2 ⇧n zwei interpolierende Polynome. Für deren Di↵erenz p := p1 p2 ist p(xi ) = 0, i = 0,..., n, d.h. das Polynom p 2 ⇧n hat n + 1 Nullstellen und verschwindet daher identisch nach Satz 2.1 (i). Abbildung 2.1: Lagrange-Polynome vom 4. Grades bei äquidistanten Stellen xj. Der Ansatz von Lagrange liefert einen wesentlich eleganteren Beweis als der naive Ansatz n X Pn (x) := c k xk k=0 in der Monom-Basis. Dieser führt auf das inhomogene lineare System n X ck xkj = yj , j = 0,..., n k=0 mit nichtsingulärer Koeffizientenmatrix (xkj )0j,kn (Übungsaufgabe !). Die Lagrange-Darstellung von Interpolationspolynomen eignet sich wegen des einfachen symmetri- schen Aufbaus gut für theoretische Zwecke. Für die praktische Berechnung ist sie weniger geeignet, da z.B. die Hinzunahme weiterer Stützstellen eine völlige Neuberechnung erfordert. Ein weiteres Pro- blem ist, dass die Lagrange-Polynome für große Werte von n stark anwachsen und oszillieren können, vgl. Beispiel 2.4. In diesem Sinne ist das Problem der Lagrange-Interpolation schlecht konditioniert. Beispiel 2.4 Wir erzeugen ein Polynom 9. Grades das an den Stellen datax := [1, 2, 3, 4, 5, 6, 7, 8, 9, 10] die zufällig  16 KAPITEL 2. INTERPOLATION Abbildung 2.2: Interpolation eines Polynoms 9. Grades mit Zufallsdaten erzeugten Werte datay := [1, 8, 3, 3, 5, 3, 9, 7, 7, 9] annimmt. Man erhält das in Abb. 2.2 gezeigte Po- lynom und sieht, dass das Polynom an den Intervallenden zum ”Überschwingen” neigt. Beispielsweise erhält man für x = 9.5 einen Wert P (9.5) > 25 und für P (0) = 836 (aufpassen, dies ist eine Extra- polation!), siehe Code-Schnipsel 2.1. Insofern ist dies ein Beispiel zur “Warnung”. Code-Schnipsel 2.1 (Zu Beispiel 2.4) 1 from numpy import linspace 2 import matplotlib. pyplot as plt 3 from scipy. interpolate import lagrange 4 # plotting : 5 def pl ot _l ag range_poly ( xdata , ydata , Pn ): 6 xplot = linspace (1 ,10 ,100) 7 yplot = Pn ( xplot ) 8 plt. plot ( xdata , ydata , 'x ') 9 plt. plot ( xplot , yplot , ' - ') 10 # interpolation : 11 xdata = [1 , 2 , 3 , 4 , 5 , 6 , 7 , 8 , 9 ,10] 12 ydata = [1 , 8 , 3 , 3 , 5 , 3 , 9 , 7 , 7 , 9] 13 Pn = lagrange ( xdata , ydata ) 14 pl ot _l ag ra ng e_poly ( xdata , ydata , Pn ) 15 plt. savefig ( " snippets / lagrange_random. pdf " ) 16 # point evaluations : 17 print ( " $P_n (9.5)={} , P_n (0)={} $ ". format ( Pn (9.5) , Pn (0))) Output: Pn (9.5) = 25.197555543067892, Pn (0) = 836.0 2.1.2 Newtonsche Interpolation Wir suchen jetzt nach einer für die praktische Berechnung besser geeignete Basis-Darstellung des Interpolationspolynoms. Die Idee der Newton-Darstellung besteht darin, die Funktionen kY1 nk (x) := (x xi ) 2 ⇧k , k = 0,... , n (2.6) i=0 als Basis des ⇧n zu nutzen. Das führt über den Ansatz 2.1. POLYNOMIALE INTERPOLATION 17 Abbildung 2.3: Newton-Basis {nk }k=0,1,2,3 bis zum Grad 3 bei äquidistanten Stellen xj = j mit j = 1,... , 3. n X Nn (x) := ck nk (x), k=0 auf das lineare Sytem n X ck nk (xj ) = yj , j = 0,..., n. (2.7) k=0 Wegen nk (xj ) = 0 für k > j ist die Koeffizientenmatrix eine untere Dreiecksmatrix, die ferner nicht- verschwindende Hauptdiagonalelemente hat. Damit können die gesuchten Werte c0 ,..., cn sukzessiv berechnet werden. Beispiel 2.5 Gegeben seien die Paare (0, 1), (1, 3), (3, 2). Dann ist die Newton-Basis n0 = 1, n1 = x, n2 = x · (x 1) und wir erhalten das Gleichungssystem 0 1 0 1 0 1 n0 (x0 ) n1 (x0 ) n2 (x0 ) c0 y0 @ n0 (x1 ) n1 (x1 ) n2 (x1 ) A · @ c1 A = @ y1 A n0 (x2 ) n1 (x2 ) n2 (x2 ) c2 y2 0 1 0 1 0 1 0 1 0 1 1 0 0 c0 1 c0 1 @ 1 1 0 A·@ c1 A = @ 3 A =) @ c1 A = @ 2 A 5 1 3 6 c2 2 c2 6 5 17 5 2 Damit ist das Interpolationspolynom N2 (x) = 1 · 1 + 2 · x + ( 6 )x · (x 1)= 1 + 6 x 6x. Eine alternative Vorgehensweise erhält man über die dividierten Di↵erenzen. 18 KAPITEL 2. INTERPOLATION Definition 2.6 Seien (xi , yi ) 2 R⇥R, i = 0,..., n gegeben mit xi 6= xk , i 6= k. Dann sind die dividierten Di↵erenzen Dik der Ordnung k an der Stelle xi rekursiv definiert durch Di0 := yi , i = 0,..., n, (2.8) k 1 Di+1 Dik 1 Dik := , i = 0,..., n k, k = 1,..., n. (2.9) xi+k xi Sinnvoll ist die Anordnung der dividierten Di↵erenzen in folgendem Di↵erenzenschema (hier n = 3) mit einem Aufwand von O(n2 ) wesentlichen Rechenoperationen: x0 y0 = D00 D01 x1 y1 = D10 D02 D11 D03 (2.10) x 2 y2 = D20 D12 D21 x 3 y3 = D30 Schema der dividierten Di↵erenzen im Fall n = 3 Beispiel 2.7 Gegeben seien wieder, wie in Beispiel 2.5, die Paare (0, 1), (1, 3), (3, 2). Dann lautet das Di↵erenzen- schema 0 1 2 5 1 3 6 1 2 3 2 Lemma 2.8 Für die dividierten Di↵erenzen gilt 0 1 i+k X i+k Y B 1 C Dik = yj B @ C, i = 0,..., n k, k = 0,..., n. (2.11) xj xr A j=i r=i r 6= j Man beachte, dass das ”Einzugsgebiet“ von Dik durch i,.., i + k begrenzt ist. Beweis. Wir führen den Induktionsbeweis nach der Ordnung k durch. Die Behauptung ist per Defi- nition richtig für k = 0. Sei nun die Aussage für k 1 bewiesen. Zur Abkürzung in der folgenden Umformung führen die folgende Notation ein: i+k Y i+k X 1 Qk,j i := , sodass Dik = yj Qk,j i zu zeigen ist. xj xr r=i j=i r 6= j 2.1. POLYNOMIALE INTERPOLATION 19 Über die Definition dividierter Di↵erenzen und die Induktionsannahme ergibt sich 1 k 1 Dik = (Di+1 Dik 1 ) xi+k xi 0 1 i+k X i+k X1 1 = @ yj Qki+11,j yj Qki 1,j A xi+k xi j=i+1 j=i 0 1 i+k X1 1 = @ yj (Qki+11,j Qki 1,j ) + yi+k Qki+11,i+k yi Qki 1,i A xi+k xi j=i+1 Nun sieht man leicht, dass Qki+11,j = (xj xi )Qk,j i für j 2 {i+1,.., i+k} und Qki 1,j = (xj xi+k )Qk,j i für j 2 {i,.., i+k 1} und damit im besonderen Qki+11,j Qki 1,j = ((xj xi ) (xj xi+k )) Qk,j i = (xi+k xi )Qk,j i für j 2 {i + 1,.., i + k 1}, Qki+11,i+k = (xi+k xi )Qk,i+k i und Qki 1,i = (xi xi+k )Qk,i i = (xi+k xi )Qk,i i gilt. Damit erhalten wir 0 1 i+k X1 i+k X Dik = @ yj Qk,j i A + yi Qk,i i + yi+k Qk,i+k i = yj Qk,j i j=i+1 j=i Satz 2.9 Das nach Satz 2.3 eindeutig bestimmte Interpolationspolynom lautet in Newton-Darstellung n n "k 1 # X X Y Nn (x) = D0k nk (x) ⌘ D0k (x xr ). (2.12) k=0 k=0 r=0 Beweis. Die rechte Seite von (2.12) sei mit Ñn (x) bezeichnet. Wir zeigen Nn (x) = Ñn (x) per Induk- tionsbeweis nach der Polynomordnung n und verwenden die Eindeutigkeit der Lagrange-Darstellung nach Satz 2.3, d.h. es gilt auch Pn ⌘ Nn. O↵enbar ist die Aussage für n = 0 richtig. Sei nun die Behauptung richtig für n 1 mit n 1. Dann betrachten wir die Di↵erenz dn := Pn Ñn. Wegen des rekursiven Aufbaus des Newton-Polynoms gilt n Y1 dn (x) = Pn (x) Ñn 1 (x) D0n (x xr ). r=0 Nach Q Lemma 2.8 und Satz 2.3 ergibt sich, dass die Koeffizienten vor xn in Pn und dem Term D0n nr=01 (x xr ) übereinstimmen. Daher verschwindet der Koeffizient vor xn im Polynom dn , d.h. dn 2 ⇧n 1. Nach Induktionsannahme ist Ñn 1 (xi ) = yi = Pn (xi ), i = 0,..., n 1, 20 KAPITEL 2. INTERPOLATION Qn 1 daher gilt mit D0n r=0 (xi xr ) = 0 für i = 0,..., n 1 auch, dass dn (xi ) = 0, i = 0,..., n 1. Damit muss dn 2 ⇧n 1 nach Satz 2.1 identisch verschwinden, also gilt Nn ⌘ Pn ⌘ Ñn. Beispiel 2.10 Das Newtonsche Interpolationspolynom zum Beispiel 2.7 lässt sich somit direkt aus dem Schema der dividierten Di↵erenzen ablesen, vgl. auch Beispiel 2.5: 5 N2 (x) = 1 + 2x x(x 1). ⇤ 6 Die Newton-Darstellung des Interpolationspolynoms hat gegenüber der Darstellung von Lagrange den wesentlichen Vorteil, dass bei Hinzunahme eines neuen Punktepaares (xn+1 , yn+1 ) die vorherige Rechnung nicht überflüssig wird. Wir betrachten nun den Aufwand der Funktionswertberechnung mittels des Newtonschen Interpola- tionspolynoms. Setzt man wie beim Horner-Schema in geeigneter Weise Klammern Nn (x) = cn (x x0 )(x x1 ) · · · (x xn 1) +.... + c1 (x x0 ) + c 0 = (...(cn (x xn 1) + cn 1 )(x xn 2) +... + c1 )(x x0 ) + c 0 , so erhält man einen Funktionswert durch n Multplikationen und 2n Additionen. Die Berechnung der Koeffizienten c0 ,..., cn erfordert selbst 12 n(n + 1) Divisionen und n(n + 1) Additionen. Man kann nun die Berechnung eines einzelnen Wertes des Interpolationspolynoms sogar ohne expli- zite Ermittlung der Polynomkoeffizienten vornehmen. Das entstehende Schema von Neville-Aitken ist dem der Berechnung der dividierten Di↵erenzen sehr ähnlich. Der entsprechende Algorithmus erfordert n(n + 1) wesentliche Rechenoperationen. Satz 2.11 Seien n + 1 Paare (xi , yi ), i = 0,..., n mit xi 6= xk , i 6= k und yi 2 R gegeben. Dann gilt für die eindeutig bestimmten Interpolationspolynome Pik 2 ⇧k , i = 0,..., n k, k = 0,..., n mit der Eigenschaft Pik (xj ) = yj , j = i,..., i + k die Rekursionsbeziehung (Schema von Neville-Aitken) Pi0 (x) = yi , (2.13) k 1 (x xi )Pi+1 (x) (x xi+k )Pik 1 (x) Pik (x) = , k = 1,..., n. (2.14) xi+k xi Der gesuchte Interpolationswert an der Stelle x ist pn (x) = P0n (x). Beweis. Wir führen den Induktionsbeweis nach der Polynomordnung k durch. Für k = 1 ist die Aussage richtig. Die Induktionsbehauptung sei nun bewiesen für k 1 mit k 2. Das durch die 2.1. POLYNOMIALE INTERPOLATION 21 zweite Rekursionsbeziehung (2.14) im Satz definierte Polynom pk gehört zu ⇧k. Weiter gilt nach Induktionsannahme (xj xi )yj (xj xi+k )yj pk (xj ) = = yj , j = i + 1,..., i + k 1 xi+k xi und pk (xi ) = yi , pk (xi+k ) = yi+k. Daraus folgt die Behauptung. Der gesuchte Interpolationswert ist pn (x) = P0n (x). Abbildung 2.4: Illustration des Di↵erenzenschemas von Neville-Aitken für Beispiel 2.7. Ausge- hend von den Graphen zu konstanten Polynomen in der linken Spalte (entsprechend (2.13)), enspricht jeder weitere Graph der Linearkombination zweier Polynome entsprechend des Di↵e- renzenschemas (2.14). Beispiel 2.12 Wir betrachten erneut die Wertepaare aus Beispiel 2.7. Diesesmal soll das Interpolationspolynom an der Stell x = 2 ausgewertet werden (ohne es vorab zu bestimmen). Dafür verwenden wir das Schema von Neville-Aitken 0 1 (2 0)·3 (2 1)·1 1 0 =5 (2 0)· 52 (2 3)·5 10 1 3 3 0 = 3 (2 1)·2 (2 3)·3 5 3 1 = 2 3 2 22 KAPITEL 2. INTERPOLATION 2.1.3 Interpolationsfehlerabschätzungen Wir betrachten die Interpolation einer gegebenen ( komplizierten”’) stetigen Funktion f : [a, b] ! R ” durch das ( einfachere”) Interpolationspolynom (Ln f )(·) 2 ⇧n mit der Interpolationseigenschaft ” (Ln f )(xi ) = f (xi ), i = 0,..., n. (2.15) Der Operator Ln : C[a, b] ! ⇧n mit f 7! Ln f heißt auch Lagrange-Interpolationsoperator. Verschiedentlich (z.B. bei der numerischen Integration) benötigt man Aussagen über den Interpola- tionsfehler f Ln f. Satz 2.13 Sei f 2 C n+1 [a, b]. Dann hat das Restglied Rn f := f Ln f bei der Polynominterpolation an n + 1 paarweise verschiedenen Stützstellen x0 ,..., xn 2 [a, b] die Lagrangesche Darstellung n f (n+1) (⇠) Y (Rn f )(x) = (x xk ), x 2 [a, b] (2.16) (n + 1)! k=0 mit einer von der Entwicklungsstelle x abhängigen Stelle ⇠ 2 [a, b]. Beweis. Die Aussage des Satzes ist erfüllt, wenn x selbst Stützstelle ist. Sei nun n Y !n+1 (x) := (x xk ) (2.17) k=0 und für festes, aber beliebiges x 2 [a, b], x 6= xk , k = 0,..., n eine Hilfsfunktion g : [a, b] ! R definiert durch f (x) (Ln f )(x) g(y) := f (y) (Ln f )(y) !n+1 (y) , y 2 [a, b]. (2.18) !n+1 (x) Nach den Voraussetzungen an f ist g 2 C n+1 [a, b]. Ferner hat g die n + 2 Nullstellen x, x0 ,..., xn. Der Satz von Rolle besagt nun, dass die Ableitung g 0 zwischen zwei Nullstellen von g wieder eine Nullstelle besitzt. Damit hat g 0 also n + 1 paarweise verschiedene Nullstellen auf [a, b]. Sukzessive Wiederholung dieses Arguments ergibt, dass g (n+1) eine Nullstelle ⇠ in [a, b] hat. Man berechnet (Rn f )(x) 0 = g (n+1) (⇠) = f (n+1) (⇠) (n + 1)! , !n+1 (x) und erhält die Behauptung. Da die Zwischenstelle ⇠ im allgemeinen Fall nicht explizit berechnet werden kann, schätzt man oft den Interpolationsfehler ab durch Korollar 2.14 Unter den Voraussetzungen von Satz 2.13 ist 1 kf L n f k1  k!n+1 k1 kf (n+1) k1 (2.19) (n + 1)! Q n mit !n+1 (x) := (x xi ) und kgk1 := max |g(x)|. i=0 x2[a,b] 2.1. POLYNOMIALE INTERPOLATION 23 Beispiel 2.15 Sei f (x) = sin(x) auf dem Intervall [0, 1]. Mit python berechnen wir die Interpolationspolynome L3 f (x) und L6 f (x) zu den äquidistanten Stützstellen xi = i/3, i = 0,..., 3 bzw. xi = i/6, i = 0,..., 6. Die Abbildung 2.5 zeigt jeweils die Lösungen L3 f bzw. L6 f sowie den Interpolationsfehler f L3 f bzw. f L6 f. O↵enbar ist f 2 C 1 [0, 1] mit kf (n) k1  1. Eine sehr grobe Abschätzung ist Abbildung 2.5: Interpolation Ln f (oben) und Interpolationsfehler f Ln f (unten) für f (x) = sin(x) und n = 3 (links) sowie n = 6 (rechts). k!n+1 k1  1. Damit folgt nach Satz 2.13 1 |f (x) (Ln f )(x)| . (n + 1)! Abb. 2.5 zeigt sogar kleinere Fehlerwerte. Mit wachsender Stützstellenzahl n wird er sogar ”expo- nentiell klein”. ⇤ Bemerkung 2.16 Eine grobe Abschätzung in Korollar 2.14 ergibt k!n+1 k1  (b a)n+1 , mit h := b a gilt dann hn+1 kf Ln f k1  kf n+1 k1. (n + 1) ! Korollar 2.14 verdeutlicht in Verbindung mit dem Beispiel, dass bei hinreichend oft di↵erenzierbaren (”glatten”) Funktionen die Interpolation mit Polynomen hoher Ordnung sinnvoll sein kann. Eventuell zerlegt man dazu das Intervall [a, b] in kleinere Intervalle und interpoliert f bei stetiger Verheftung an den Intervallenden stückweise durch Lagrange-Polynome. Das ist auch die Grundidee der in An- wendungen sehr wichtigen Finite-Elemente-Methoden höherer Ordnung. 24 KAPITEL 2. INTERPOLATION Zu beachten ist jedoch, dass auch die Größen kf (n+1) k1 nicht beliebig stark bezüglich n anwach- sen sollten. Runge untersuchte das Beispiel der auf [a, b] = [ 1, 1] beliebig oft di↵erenzierbaren Funktion 1 f (x) =. 1 + 25x2 Man beobachtet bereits für moderate Werte von n die Divergenz der Interpolationspolynome an den Intervallenden (vgl. Übungsaufgabe). Bemerkung 2.17 Korollar 2.14 legt nahe, bei Wahlmöglichkeit der Stützstellen diese so festzulegen, dass n Y max (x xk ) ! M in.! (2.20) x2[a,b] k=0 In der Approximationstheorie wird gezeigt, dass dabei die Stützstellen ✓ ◆ a+b b a 2(n i) + 1 xi = + cos ⇡ , i = 0,..., n (2.21) 2 2 2(n + 1) entstehen (Nullstellen der Tschebychev-Polynome). Bei dieser Wahl erreicht man auch eine wesentlich bessere Approximation für das Beispiel von Runge (vgl. Übungsaufgabe). Bei der Interpolation haben wir bisher nur die Qualität der Auswertung im Intervall betrachtet. Tatsächlich kann das Interpolationspolynom aber auch zum Auswerten von Ableitungen herangezogen werden. Jedoch gilt, dass für jede Ableitung “eine Ordnung verloren geht”, d.h. dass sich die h-Potenz (h: Intervallgröße) um eins reduziert, wie wir im folgenden Lemma sehen: Lemma 2.17b (Abschätzungen für die Approximation der Ableitungen1 ) Sei f 2 C n+1 [a, b]. Dann gilt für die m-te Ableitung des Restglieds Rn f := f Ln f , m 2 {0,.., n + 1}, bei der Polynominterpolation an n + 1 paarweise verschiedenen Stützstellen x0 ,..., xn 2 [a, b] die Abschätzung hn+1 m k(f Ln f )(m) k1  kf (n+1) k1 (n + 1 m)! Beweis. (s. Übung.) 2.1.4 Konvergenz von Interpolationspolynomen Wir betrachten nun die Konvergenz von Interpolationspolynomen bei wachsender Stützstellenzahl. O↵enbar gilt Satz 2.18 Seien f 2 C 1 [a, b] und kf (n) k1  M für alle n = 0, 1,.... Dann konvergiert der Interpolationsfehler Rn f für n ! 1 gleichmäßig auf [a, b] gegen Null. n+1 Beweis. Die Behauptung folgt aus |(Rn f )(x)|  M (b(n+1)! a) ! 0, n ! 1. 1 Der Index 2.17b wurde hier nur reingeschummelt, um die nachfolgende Zählung nicht zu verändern. 2.1. POLYNOMIALE INTERPOLATION 25 Dieser Satz ist zwar für viele Standardfunktionen, jedoch in der Regel nicht auf praktisch interessieren- de Funktionen, anwendbar. Insbesondere ist die starke Glattheitsforderung an die zu interpolierende Funktion stark einschränkend. Das folgende Beispiel zeigt, dass im Fall lediglich stetiger Funktionen der Fehler bei wachsender Stützstellenzahl sogar divergieren kann. Beispiel 2.19 Sei ⇢ x sin ⇡x , x 2 (0, 1], f (x) := 0, x = 0. 1 Mit xk = k+1 , k = 0,..., n ist wegen f (xk ) = 0, k = 0,..., n o↵enbar (Ln f )(x) = 0 für alle n 2 N. Die Folge der Interpolationspolynome konvergiert bei dieser Stützstellenwahl nur an den Stellen xk , k 2 N0 gegen die Funktion f. Zur Verdeutlichung führen wir noch ohne Beweis folgende Sätze über die Konvergenz von Interpo- lationspolynomen an. Satz 2.20 (Marcinkiewicz) (n) Für jede Funktion f 2 C[a, b] gibt es eine Folge von Stützstellen (xk ), k = 0,..., n und n 2 N0 der- (n) art, dass die entsprechende Folge (Ln f ) von Interpolationspolynomen Ln f 2 ⇧n mit (Ln f )(xk ) = (n) f (xk ), k = 0,..., n auf [a, b] gleichmäßig gegen f konvergiert. Satz 2.21 (Faber) (n) Für jede Folge (xk ) existiert eine Funktion f 2 C[a, b] so, dass die zugehörige Folge (Ln f ) von Interpolationspolynomen auf [a, b] nicht gleichmäßig gegen f konvergiert. Bei geeigneter Wahl der Stützstellen gilt jedoch wenigstens Satz 2.22 Sei f 2 C 1 [a, b]. Dann konvergiert die Folge der mit den Tschebychev-Stützstellen (vgl. Bemer- kung 2.17) gebildeten Interpolationspolynome gleichmäßig auf [a, b] gegen f. Insgesamt ergibt sich aus den Überlegungen dieses Abschnitts, dass die Interpolation mit Polyno- men hohen Grades im allgemeinen Fall (insbesondere bei geringer Glattheit der zu interpolierenden Funktionen) nicht sinnvoll ist. Wir werden im Kapitel 2.3 über die Spline-Interpolation sehen, wie durch stückweise polynomiale Interpolation unter relativ geringen Glattheitsforderungen an die zu interpolierende Funktion gleichmäßige Konvergenz erzielt werden kann. 2.1.5 Verallgemeinerung Wir betrachten die Interpolation mittels allgemeinerer Funktionensysteme. Definition 2.23 Ein m dimensionaler Unterraum U ⇢ C[a, b] heißt unisolvent bezüglich der paarweise verschiedenen Stützstellen x1 ,..., xm 2 [a, b], wenn jede Funktion u 2 U mit den Nullstellen u(xi ) = 0, i = 1,...m identisch verschwindet. 26 KAPITEL 2. INTERPOLATION Im vorliegenden Kapitel hatten wir den Fall U = ⇧n [a, b] ⇢ C[a, b] mit verschiedenen Basisfunktionen (Monom-Basis, Basis mit Lagrange- bzw. Newton-Polynomen) und m = n + 1 besprochen. Satz 2.24 Seien der m dimensionale Unterraum U ⇢ C[a, b] bezüglich der paarweise verschiedenen Stützstellen x1 ,..., xm 2 [a, b] unisolvent und m Werte y1 ,..., ym 2 R gegeben. Dann existiert genau eine Funk- tion u 2 U mit der Interpolationseigenschaft u(xi ) = yi , i = 1,..., m. Beweis. Sei U = span{ 1 ,..., m }. Mit dem Ansatz m X u= ck k k=1 führt die Interpolationsaufgabe auf das lineare Gleichungssystem m X ck k (xi ) = yi , i = 1,..., m. (2.22) k=1 Wegen der Bedingung der Unisolvenz an U hat das zugehörige homogene System nur die triviale Lösung. Folglich ist das Gleichungssystem eindeutig lösbar. 2.2 Trigonometrische Interpolation Bei Anwendungen in der Signal- und Steuerungstechnik treten in der Regel periodische Daten bzw. periodische Funktionen mit der Eigenschaft 9T > 0 : f (x + T ) = f (x), 8x 2 R (2.23) auf. Unter dem Aspekt der modernen Kommunikationstechnik kommt es hier vor allem auf eine sehr effiziente Auswertung derartiger Daten an. Algebraische Polynome (vgl. Kapitel 2.1) sind nicht zur Interpolation derartiger Daten geeignet, da sie nicht periodisch sind. Wir wenden uns daher trigonometrischen Polynomen zu. Bei äquidistanter Verteilung der Interpolationsstützstellen spricht man auch von diskreter Fourier-Transformation. 2.2.1 Trigonometrische Polynome Wir werden ohne Beschränkung der Allgemeinheit annehmen, dass für die Periodenlänge gilt T = 2⇡. Dann ist die Betrachtung der folgenden Funktionensysteme sinnvoll. Definition 2.25 Die Elemente des Raumes Tn := span{1, cos x, sin x,..., cos(nx), sin(nx)} (2.24) heißen trigonometrische oder Fourier-Polynome vom Grad kleiner oder gleich n 2 N0. 2.2. TRIGONOMETRISCHE INTERPOLATION 27 Die Begri↵sbildung wird durch folgende Überlegung gerechtfertigt: Über die Additionstheoreme für trigonometrische Funktionen folgt, dass mit p1 2 Tn1 und p2 2 Tn2 für das Produkt gilt p1 p2 2 Tn1 +n2. (Übungsaufgabe !) Nachfolgend wollen wir eine funktionalanalytische Charakterisierung des Raumes Tn vornehmen und einen Zusamenhang zwischen Fourier- und algebraischen Polynomen herstellen. Dazu beweisen wir zunächst das folgende Resultat. Lemma 2.26 Es gelten folgende Aussagen: (i) Hat ein trigonometrisches Polynom aus Tn mehr als 2n paarweise verschiedene Nullstellen im Periodizitätsintervall [0, 2⇡), so verschwindet es dort identisch. (ii) Die Funktionen cos(kx), k = 0,... , n und sin(kx), k = 1,... , n sind linear unabhängig auf dem Raum C[0, 2⇡]. Beweis. (i) Wir betrachten ein trigonometrisches Polynom pn 2 Tn der Form X n 1 pn (x) = a0 + [ak cos(kx) + bk sin(kx)] (2.25) 2 k=1 und gehen nun zu einer Formulierung in der Menge C über. Dabei ist i die imaginäre Einheit. Unter Verwendung von 1 1 ck := (ak ibk ), c k := (ak + ibk ), k = 0,..., n 2 2 und mit b0 = 0 (und damit c0 = 12 a0 ) gilt ak = ck + c k, bk = i(ck c k ), k = (0, )1..., n. Aus der Eulerschen Formel eit = cos t + i sin t ergibt sich die zu (2.25) äquivalente komplexe Formulierung n h X i n X pn (x) = c0 + (ck + c k ) cos(kx) + i(ck c k ) sin(kx) = ck eikx. (2.26) | {z } k=1 k= n =ick sin(kx)+ic k sin( kx) Mit der Substitution z = eix und der Bezeichnung n X q2n (z) = ck z n+k k= n erhalten wir n pn (x) = z q2n (z). Wir nehmen an, dass das trigonometrische Polynom pn mehr als 2n paarweise verschiedene Nullstellen in [0, 2⇡) hat. Dann hat das algebraische Polynom q2n 2 ⇧2n mehr als 2n paarweise 28 KAPITEL 2. INTERPOLATION verschiedene Nullstellen auf dem Einheitskreis in C, denn die Funktion t 7! eit bildet das Intervall [0, 2⇡) bijektiv auf den Einheitskreis in C ab. Nach Satz (2.1) muss das algebraische Polynom q2n identisch verschwinden. Wegen (2.26) muss auch das trigonometrische Polynom pn identisch verschwinden. (ii) Wir nehmen an, dass X n 1 a0 + [ak cos(kx) + bk sin(kx)] = 0, x 2 [0, 2⇡]. 2 k=1 Das trigonometrische Polynom mit den Koe↵zienten a0 , a1 , b1 ,... , an , bn hat o↵enbar mehr als 2n paarweise verschiedene Nullstellen in [0, 2⇡). Nach Aussage (i) müssen dann aber die Koeffizienten sämtlich verschwinden. Als Folgerung aus Lemma 2.2 erhalten wir folgendes wesentliche Ergebnis. Satz 2.27 Der lineare Raum Tn der trigonometrischen Polynome vom Grad kleiner oder gleich n hat die Di- mension 2n + 1 und ist unisolvent bezüglich 2n + 1 paarweise verschiedener Stützstellen aus dem Intervall [0, 2⇡). 2.2.2 Trigonometrische Interpolation Definition 2.28 Seien 2n + 1 paarweise verschiedene Stützstellen x0 ,..., x2n im Intervall [0, 2⇡) sowie 2n + 1 Werte y0 ,..., y2n 2 R vorgegeben. Bei der trigonometrischen Interpolation sucht man ein trigonometrisches Polynom pn 2 Tn mit der Interpolationseigenschaft pn (xj ) = yj , j = 0,..., 2n. Nachfolgender Satz zeigt, dass diese Aufgabe wohlgestellt ist. Satz 2.29 Die Aufgabe der trigonometrischen Interpolation ist eindeutig lösbar. Mit Hilfe der Lagrange-Polynome 2n Y sin 12 (x xj ) lk (x) = , k = 0,..., 2n (2.27) j =0 sin 12 (xk xj ) j 6= k lautet die Lagrange-Darstellung des Interpolationspolynoms 2n X pn (x) = yk lk (x). (2.28) k=0 Beweis. 2.2. TRIGONOMETRISCHE INTERPOLATION 29 (i) Zunächst stellen wir die Existenz und Eindeutigkeit der Interpolationsfunktion fest, da wegen Satz 2.27 der Satz 2.24 anwendbar ist. (ii) Die konkrete Gestalt von pn folgt wegen lk (xj ) = kj und lk 2 Tn , k = 0,..., 2n. Die Aussage lk 2 Tn ergibt sich dabei aus dem Additionstheorem ✓ ◆ x x0 x x1 1 x1 x0 1 x1 + x0 sin sin = cos cos x , 2 2 2 2 2 2 d.h. lk ist ein Produkt aus n trigonometrischen Polynomen vom Grad 1. Abbildung 2.6: Trigonometrische Lagrange-Polynome bei äquidistanten Stellen und n = 0, 1, 2. Wir hatten bereits im Fall der Interpolation mit algebraischen Polynomen gesehen, dass die Lagrange- Basis für praktische Zwecke wenig geeignet ist. Dies gilt natürlich auch im Fall der trigonometrischen Interpolation. Daher gehen wir künftig von dem Ansatz (2.25) bzw. (2.26) aus. Das Problem besteht somit in der effizienten Bestimmung der Koeffizienten ak , bk bzw. ck. Zur Vereinfachung der Darstellung betrachen wir nun den wichtigen Spezialfall äquidistanter Stützstellen. Sei zunächst eine ungerade Anzahl von Stützstellen N = 2n + 1 angenommen mit 2⇡ xj = j , j = 0,..., 2n. (2.29) N Dann gilt die Summenformel 2n X ⇢ ikxj N, k = 0, e = , (2.30) 0, k = ±1,..., ±2n j=0 denn für eixk 6= 1 ergeben die geometrische Summe und die Eulersche Formel 2n X 2n X 2n X j 1 eiN xk eikxj = eijxk = eixk = = 0. 1 eixk j=0 j=0 j=0 Zur Vereinfachung der Berechnung setzen wir das trigonometrische Interpolationspolynom in kom- plexer Form an X n pn (x) = ck eikx. (2.26) k= n 30 KAPITEL 2. INTERPOLATION Aus der Interpolationsforderung pn (xj ) = yj , j = 0,..., 2n schließen wir auf das eindeutig lösbare lineare Gleichungssystem n X ck eikxj = yj , j = 0,..., 2n. k= n Zur Lösung des Systems multiplizieren wir mit e irxj , r = n,... , n, summieren über j, vertauschen die Summationsfolge und benutzen Beziehung (2.30) 2n X n X 2n X irxj yj e = ck ei(k r)xj = N cr j=0 k= n j=0 und daher (mit 2n = N 1) N 1 1 X ikxj ck = yj e , k= n,..., n. (2.31) N j=0 Die Formeln (2.25) und (2.26) zeigen den Zusammenhang zwischen reeller und komplexer Formulie- rung des trigonometrischen Polynoms, insbesondere ist ak = ck + c k, bk = i(ck c k ), k = 0,... , n. Daraus ergibt sich folgendes Resultat. Satz 2.30 Sei N = 2n + 1. Es existiert genau ein trigonometrisches Polynom X n 1 pn (x) = a0 + [ak cos(kx) + bk sin(kx)] 2 k=1 mit der Interpolationseigenschaft ✓ ◆ 2⇡ pn j = yj , j = 0,..., 2n. N Für die Koeffizienten gelten die Formeln N 1 ✓ ◆ 2 X 2⇡ ak = yj cos jk , k = 0,..., n (2.32a) N N j=0 N X1 ✓ ◆ 2 2⇡ bk = yj sin jk , k = 1,..., n. (2.32b) N N j=0 Im Falle einer geraden Anzahl der Stützstellen N = 2n gibt es im äquidistanten Fall mit 2⇡ xj = j. j = 0,..., 2n 1 (2.33) N zu den N +1 = 2n+1 Freiheitsgraden eines trigonometrischen Polynoms aus dem Raum Tn zunächst nur N = 2n Interpolationsbedingungen. Man kann jedoch die Basisfunktion sin nx weglassen, da sie in den Stützstellen nur Nullstellen hat. Dann gilt 2.2. TRIGONOMETRISCHE INTERPOLATION 31 Satz 2.31 Sei N = 2n. Es existiert genau ein trigonometrisches Polynom X n 1 1 1 pn (x) = a0 + [ak cos(kx) + bk sin(kx)] + an cos(nx) 2 2 k=1 mit der Interpolationseigenschaft ✓ ◆ 2⇡ pn j = yj , j = 0,..., 2n 1. N Die Koeffizienten sind N 1 ✓ ◆ 2 X 2⇡ ak = yj cos jk , k = 0,..., n, (2.34a) N N j=0 N 1 ✓ ◆ 2 X 2⇡ bk = yj sin jk , k = 1,..., n 1. (2.34b) N N j=0 Beweis. Der lineare Raum T̃n der trigonometrischen Polynome aus Tn mit verschwindendem Koeffizi- enten bn = 0 vor sin(nx) hat o↵enbar die Dimension 2n. Wir betrachten nun ein Element p 2 T̃n mit den Nullstellen p(xj ) = 0, j = 0,..., 2n 1. Ferner wird der Punkt x⇤ mit x⇤ 6= xj , j = 0,..., 2n 1 gewählt und sin(nx) Q(x) := p(x) p(x⇤ ) sin(nx⇤ ) gesetzt. Dann ist Q 2 Tn und hat die 2n + 1 Nullstellen xj , j = 0,...2n 1 sowie x⇤. Daher muss Q verschwinden und p ist ein Vielfaches von sin(nx). Nach Definition verschwindet dann auch p 2 T̃n. Daher ist T̃n unisolvent bezüglich der äquidistanten Stützstellen xj , j = 0,..., 2n 1. Existenz und Eindeutigkeit der Interpolationsfunktion ergeben sich dann aus Satz 2.24. Die gesuchte Darstellung gewinnt man analog zum Beweis von Satz 2.30. Man spricht im Fall äquidistanter Stützstellenverteilung von der diskreten Fourier-Transformation. Wir wollen diesen Punkt etwas genauer betrachten: Bemerkung 2.32 Man kann zeigen (vgl. Übungsaufgabe), dass eine gegebene Funktion f 2 C[0, 2⇡] in der Klasse der trigonometrischen Polynome vom maximalen Grad n bezüglich der L2 (0, 2⇡)-Norm ⇣ Z 2⇡ ⌘1 2 kgkL2 (0,2⇡) := [g(x)]2 dx 0 am besten durch das Polynom 1 Xh n i (Pn f )(x) = ã0 + ãk cos(kx) + b̃k sin(kx) 2 k=1 mit den Fourier-Koeffizienten (für k = 0,... , n) Z Z 1 2⇡ 1 2⇡ ãk = f (x) cos(kx) dx, b̃k = f (x) sin(kx) dx (2.35) ⇡ 0 ⇡ 0 32 KAPITEL 2. INTERPOLATION approximiert wird. Man erhält diese Formeln auch, wenn man in den Formeln (2.32a), (2.32b) für ak bzw. bk formal den Grenzübergang n ! 1 ausführt. Andererseits erhält man die Formeln (2.32a), (2.32b) aus den Formeln (2.35) bei Anwendung der zusammengesetzten numerischen Integration mittels der sogenannten Trapezregel, die wir im Kapitel 3.1 behandeln werden. Der aufgezeigte Zusammenhang verdeutlicht die Begri↵sbildung ”diskrete Fourier-Transformation”. 2.2.3 Berechnung der Fourier-Koeffizienten Wir betrachten hier nur die tragende Idee zur Berechnung der komplexwertigen Fourier-Koeffizienten. Seien dazu yj die gegebenen Funktionswerte, dann suchen wir die Koeffizienten ck des trigonome- trischen Polynoms, also N 1 1 X 2⇡i kj ck = yj e N , k = 0,..., N 1. N j=0 2⇡i N = 1. Wir Sei nun !N = e N die N -te komplexe Einheitswurzel. Dann gilt o↵ensichtlich !N definieren die Matrix ⇣ ⌘N 1 kj VN = !N , k,j=0 und sehen, dass sich die Koeffizienten ck mit Hilfe von 0 1 0 1 c0 y0 B.. C 1 B. C @. A = VN @.. A , N cN 1 yN 1 berechnen lassen. Die Matrix VN ist eine sogenannte DFT-Matrix, und die Berechnung der Fourier- Koeffizienten entspricht einer diskreten Fourier-Transformation (siehe Bemerkung 2.32). Lemma 2.33 Die DFT-Matrix VN ist unitär, d.h. es gilt VN1 = 1 N VN , wobei VN die konjugiert-transponierte Matrix von VN ist. Damit sehen wir auch, dass V N V N = N IN. Beweis. Übungsbeispiel. Für die Berechnung der Fourier-Koeffizienten greifen wir also auf ein Matrix-Vektor-Produkt zurück. Die Berechnung der Fourier-Koeffizienten ist also mit O(N 2 ) Operationen verbunden. Bei großen Werten von n möchte man diesen Aufwand weiter reduzieren. Eine dafür geeignete Methode ist die schnelle Fourier-Transformation (FFT - fast Fourier transformation) nach Cooley/ Tukey (1965). Eine geschickte Ausnutzung von Symmetrien der Einheitswurzeln in C im Falle von n = 2s ist dabei der Schlüssel. Eine gute Gesamtdarstellung findet man z.B. im Lehrbuch , Kapitel 53 oder in dem Artikel von H.R. Schwarz: Elementare Darstellung der schnellen Fourier-Transformation in Computing 18 (1977), 107-116. Im Folgenden betrachten wir den Fall einer geraden Anzahl von Stützstellen N = 2n.  2.2. TRIGONOMETRISCHE INTERPOLATION 33 Satz 2.34 (Schnelle-Fourier-Transformation) 2⇡i Sei N = 2n mit n 2 N, und sei !N = e( N ) sowie ⇠ = (!N 2 ). Für die Berechnung der (unnormier- ten) Fourier-Koeffizienten N X1 kj ck := yj ! N , k = 0,..., N 1, j=0 gelten die Beziehungen N/2 1 X c2l = gj ⇠ jl , mit gj = (yj + yj+N/2 ), l = 0,..., N/2 1 j=0 N/2 1 X j c2l+1 = hj ⇠ jl , mit hj = (yj yj+N/2 )!N , l = 0,..., N/2 1. j=0 Beweis. Per Definition der Koeffizienten folgt N N/2 1 X1 2lj X 2lj 2l(j+N/2) c2l = yj !N = yj ! N + yj+N/2 !N j=0 j=0 N/2 1 N/2 1 X 2lj 2lj N l X 2 jl = yj !N + yj+N/2 !N !N = gj (!N ) , j=0 j=0 N = 1 verwendet haben. Der Beweis für c wobei wir im letzten Schritt die Eigenschaft !N 2l+1 erfolgt analog. Wir sehen also, dass sich, die DFT mit geradem N in zwei DFTs mit N/2 zerlegen lässt. Man beachte, dass wir bei der Berechnung der Koeffizienten von der Normalisierung durch N abgesehen haben. Eine sehr rudimentäre Implementierung der FFT in python findet sich im nachfolgenden Code- Schnipsel 2.2. Hierbei ist angenommen, dass das Datenfeld 2s , s 2 N lang ist. Code-Schnipsel 2.2 (FFT) 1 from math import pi 2 from cmath import exp 3 4 def fft ( y ): 5 N = len ( y ) 6 c = [0 for i in range ( N )] 7 g = [ None for i in range ( N //2)] 8 h = [ None for i in range ( N //2)] 9 if N == 1: 10 return y 11 else : 12 omega = exp ( -2* pi *1 j / N ) 13 for l in range ( N //2): 14 g [ l ] = y [ l ] + y [ l + N //2] 15 h [ l ] = ( y [ l ] - y [ l + N //2]) * omega ** l 16 c [::2] = fft ( g ) 17 c [1::2] = fft ( h ) 18 return c Für die Kosten zur Berechnung der Fourier-Koeffizienten ergibt sich damit. 34 KAPITEL 2. INTERPOLATION Lemma 2.35 Die Berechnung der Fourier-Koeffizienten eines trigonometrischen Polynoms vom Grad n mit N = 2n Stützstellen benötigt O(N log N ) Operationen. Beweis. Seien AN die Kosten zur Berechnung der Fourier-Koeffizienten eines trigonometrischen Polynoms vom Grad n mit N = 2n Stützstellen. Dann gilt AN = 2AN/2 + cN, (2.36) wobei cN (mit einer Konstanten c) die Kosten für die Berechnung der Koeffizienten gj und hj in Theorem 2.34 sind. Mit A1 = c erhalten wir die Lösung AN = cN (log2 N + 1) denn mit (2.36) folgt (2.36) N AN = cN (log2 N + 1) = 2(c (log2 N/2 + 1)) + cN = 2N (log2 N 1 + 1) + cN = AN. 2 Die Behauptung folgt mit O(N log(N + 1)) = O(N log N ). Beispiel 2.36 Die folgende Überlegung verdeutlicht den erheblich reduzierten Aufwand bei der FFT. Sei n = 28 = 256. Dann ist n2 = 65.536 und n · log2 n = 2.048, d.h. der Aufwand der FFT liegt um den Faktor 32 niedriger als bei der Anwendung mittels einer Matrix-Vektor-Multiplikation. 2.2.4 Konvergenz trigonometrischer Polynome Konvergenzuntersuchungen für die Folge (pn )n der trigonometrischen Interpolationspolynome an eine gegebene periodische Funktion f gestalten sich komplizierter als entsprechende Untersuchungen bei algebraischen Interpolationsproblemen (vgl. Kap. 1). Es gilt jedoch das folgende Resultat. Satz 2.37 Sei f 2 C(R) eine gegebene 2⇡ periodische Funktion. Dann konvergiert die Folge der trigonome- trischen Interpolationspolynome pn an f in der L2 -Norm gegen f , d.h. lim kpn f kL2 (0,2⇡) = 0. n!1 Ist darüber hinaus sogar f 2 C 1 (R), so konvergiert die Folge in der Maximumnorm, d.h. lim kpn f k1 = 0 n!1 bzw. f (x) = limn!1 pn (x) für alle x 2 [0, 2⇡]. Ohne Beweis. (Beim Beweis benutzt man maßgeblich den Approximationssatz von Weierstraß für trigonometrische Polynome.) Insbesondere ist das Konvergenzresultat in der L2 Norm im Zusammenhang mit der Minimaleigen- schaft von Fourier-Reihen (vgl. Bemerkung 2.32) als Grenzfall für n ! 1 zu sehen. Man kann die An- wendung der Fourier-Approximation auch auf solche Funktionen f erweitern für die kf kL2 (0,2⇡) noch endlich ist. Diese Funktionen bilden den Lebesgue-Raum L2 (0, 2⇡), der durch ”Vervollständigung” von C[0, 2⇡] bezüglich der Norm k · kL2 (0,2⇡) ensteht. Das folgende Beispiel verdeutlicht, dass die (punktweise) Konvergenz der Fourier-Polynome für n ! 1 nicht sehr schnell ist, wenn die gegebene Kurve nicht mehr stetig ist.  2.2. TRIGONOMETRISCHE INTERPOLATION 35 Beispiel 2.38 Wir zeigen als Beispiel für die ”Sägezahn”-Impulskurve f (x) = x, 0  x < 1; f (x + 1) = f (x), x2R die mit python berechneten Fourier-Ansätze mit n = 4, 16, 64. Hierfür benutzen wir die reelle FFT von numpy (rfft), siehe Code-Schnipsel 2.3. Abb. 2.7 verdeutlicht, dass sich die Approximation mit wachsendem n nur langsam der vorgegebenen Funktion nähert, insbesondere in deren Unstetigkeits- punkten. Bei den lokalen Oszillationen in deren Umgebung spricht man auch vom ”Gibbs-Phänomen”. Abbildung 2.7: Fourier-Approximation der ”Sägezahnkurve” mit n = 4, 16, 64. Code-Schnipsel 2.3 (Zu Beispiel 2.38) 1 from numpy import linspace 2 from numpy. fft import rfft , irfft 3 import matplotlib. pyplot as plt 4 x = linspace (0 ,1 ,2**10 , endpoint = False ) 5 y = x % 1 6 def s e l e c t _ f i r s t _ f o u r i e r _ m o d e s (y , N ): 7 fourier_trans = rfft ( y ) 8 M = len ( fourier_trans ) 9 fourier_trans [ N +1: M ] = 0 10 return irfft ( fourier_trans ) 11 for M in [0 ,1 ,2 ,3 ,4 ,5 ,6 ,7 ,8 ,16 ,32 ,64 ,128]: 12 yF = s e l e c t _ f i r s t _ f o u r i e r _ m o d e s (y , M ) 13 plt. clf () 14 for i in range (3): 15 plt. plot ( x +i , y , ' -b ') 16 plt. plot ( x +i , yF , ' -r ') 17 plt. savefig ( " snippets / fourier_sae ge za hn _ " + str ( M )+ ". pdf " ) Für den Fall hinreichender glatter periodischer Funktionen zitieren wir ohne Beweis noch folgendes Resultat über exponentielle Konvergenz. Satz 2.39 Sei f 2 C k (R) 2⇡ periodisch mit k 2 N und sei pn das trigonometrische Interpolationspolynom nach Satz 2.30 bzw. 2.31. Dann gibt es eine nur von k abhängige Konstante Ck mit ✓Z 2⇡ h i ◆ 12 k 2 (k) 2 kf p n kL 2  C k n |f (x)| + |f (x)| dx , 8n 2 N0. 0 Insbesondere konvergiert für f 2 C 1 (R) der Fehler kf pn kL2 für n ! 1 schneller als jede Potenz n k bei festem k 2 N gegen Null. 36 KAPITEL 2. INTERPOLATION 2.3 Spline-Interpolation Im Kapitel 2.1 hatten wir festgestellt, dass Interpolationspolynome bei Vergrößerung der Zahl n der Stützstellen, d.h. genauer im Grenzprozeß n ! 1, nicht notwendig gegen die zu interpolierende Funktion f konvergieren. Einen Ausweg bietet die stückweise polynomiale Interpolation durch Splines. Sie bilden auch eine wichtige Grundlage für die (stückweise) numerische Integration (vgl. Kapitel 3.1) sowie für die Diskretisierung von Di↵erentialgleichungen. 2.3.1 Räume von Spline-Funktionen Definition 2.40 Sei a = x0 < x1 <... < xn = b eine Unterteilung des Intervalls [a, b]. Dann heißt eine Funktion s : [a, b] ! R Spline m-ten Grades, falls (i) s 2 Cm 1 [a, b], (ii) s [xj 1 ,xj ] 2 ⇧m [xj 1 , xj ], j = 1,..., n. Die Menge der Splines m ter Ordnung zu einer gegebenen Unterteilung von [a, b] in n Teilintervalle bezeichnen wir mit Snm ([a, b]). Beispiel 2.41 (Lineare Splines (m = 1)) Splines 1. Ordnung sind (global) stetige und stückweise lineare Funktionen. Bei gegebener Funktion f 2 C[a, b] lauten die Interpolationforderungen s(xj ) = f (xj ), j = 0,..., n, d.h. die Spline-Funktion entsteht einfach durch stetige Verheftung der auf den Teilintervallen stückweise linearen Interpolationspolynome in den Punkten (xj , f (xj )). Auf jedem Teilintervall [xj 1 , xj ] gilt nach Satz 2.13 die Fehlerabschätzung 1 |(s f )(x)|  |f 00 (⇠)|h2j , ⇠ 2 [xj 1 , xj ], hj := xj xj 1. 8 Damit konvergiert die Spline-Funktion 1. Ordnung auf [a, b] zumindest für f 2 C 2 [a, b] gleichmäßig gegen f , wenn h := maxj=1,...,n hj ! 0. Bei manchen Anwendungen (z.B. Design-Problemen) sind Splines 1. Ordnung nachteilig, da ein “glatterer” Übergang an den Teilpunkten erforderlich wäre. Einen Kompromiß zwischen Glattheit und Rechenaufwand bieten kubische Splines. Beispiel 2.42 (Kubische Splines (m = 3)) Nach Definition gilt s 2 C 2 [a, b] und s 2 ⇧3 [xj 1 , xj ] für j = 1,..., n. Dem englischen Wort spline entspricht etwa das deutsche Wort Strak. Das ist ein Werkzeug aus dem Schi↵bau zur Führung glatter Kurven durch vorgegebene Punkte, d.h. der Strak wird dort durch Lager fixiert. In der Mechanik leitet man für die Balkenbiegung u(x) die Gleichung u(4) (x) = 0 ab. Zwischen den Lagerpunkten wird der Strak also durch ein Polynom 3. Grades beschrieben. Eine mathematische Darstellung der Biegelinie u = u(x) erhalten wir dann unter Beachtung der Stetigkeit der Auslenkung u, der Biegung u0 sowie des Biegemomentes u00 in den Lagerpunkten. 2.3. SPLINE-INTERPOLATION 37 In Abb. 2.8 vergleichen wir den in Abschn. 2.1.1 betrachteten Datensatz (bis auf Skalierung, (vgl. Abb. 2.2) die mit der numpy-Funktion make interp spline für m = 1, 2, 3 erzeugten Spline-Funktion mit dem globalen Interpolationspolynom 9. Ordnung. Die Spline-Funktion ist zwar an den Kopp- lungspunkten weniger glatt, vermeidet aber die Randoszillationen des globalen Polynoms. Ferner Abbildung 2.8: Interpolation von Zufallsdaten. Spline-Interpolation von 1. bis 3. Ordnunge (links) und Vergleich Spline-Interpolation vom Grad 3 mit Lagrange-Interpolation vom Grad 10 (rechts) betrachten wir in Abb. 2.8 die Interpolation durch Splines 1., 2. und 3. Ordnung. Der ”glattere” Ver- lauf für kubische Splines ist o↵ensichtlich. Vergleicht man ferner die Splines mit m = 3 und m = 9 (hier nicht gezeigt), so ist die Spline-Funktion mit m = 9 zwar ”glatter”, sie tendiert aber ebenso wie das globale Polynom 9. Grades zu Randoszillationen. Auch dies spricht für die Wahl m = 3. ⇤ Zur Ermittlung einer Basis des Raumes Snm ([a, b]) der Splines der Ordnung m definieren wir über die Funktion ⇢ m m x , x 0 x+ := 0, x < 0. das System der Kardinalsplines ⇢ k (x) := (x x0 ) k , k = 0,... , m (2.37) j (x) := (x xj ) m +, j = 1,... , n 1. Abbildung 2.9: Kardinalsplines für xi = i, i = 0,.., 4 (n=4) und m = 0, 1, 2 (von links nach rechts). Lemma 2.43 Die n + m Funktionen des Systems (2.37) sind linear unabhängig. 38 KAPITEL 2. INTERPOLATION Beweis. Sei m X n X1 ak (x x0 ) k + bj (x xj ) m + = 0, x 2 [a, b]. k=0 j=1 Nach Definition ist m X ak (x x0 )k = 0, x 2 [x0 , x1 ], k=0 damit ak = 0, k = 0,... , m. Weiter gilt b1 (x x1 )m + = 0, x 2 [x1 , x2 ] und daher b1 = 0. Sukzessiv schließen wir auf bj = 0, j = 1,..., n 1. Daraus folgt die lineare Unabhängigkeit der Funktionen des Systems (2.37). Der nachfolgende Satz zeigt, dass die Funktionen des Systems (2.37) eine Basis des Snm ([a, b]) bilden. Jede Funktion s 2 Snm ([a, b]) hat dann die Kardinalspline-Darstellung m X n X1 k s(x) = ak (x x0 ) + bj (x xj ) m +, x 2 [a, b]. (2.38) k=0 j=1 Satz 2.44 Der Raum Snm ([a, b]) der Splines ist ein linearer Raum der Dimension m + n. Die Funktionen des Systems (2.37) bilden eine Basis von Snm ([a, b]). Beweis. 1. Wir zeigen die Gültigkeit der Darstellung (2.38) durch vollständige Induktion über die Teilin- tervalle, d.h. m X i 1 X k s(x) = ak (x x0 ) + bj (x xj ) m +, x 2 [x0 , xi ], i = 1,..., n 1. k=0 j=1 Der Induktionsanfang für i = 1 ist o↵enbar richtig wegen s 2 ⇧m [x0 , x1 ]. Sei nun die Gültigkeit der Darstellung bewiesen für eine Zahl i 2 {1,... , n 1}. Dann verschwindet die Di↵erenz m X i 1 X d(x) := s(x) ak (x x0 ) k bj (x xj ) m + k=0 j=1 auf dem Intervall [x0 , xi ]. Ferner ist d 2 ⇧m [xi , xi+1 ]. Weiter gilt dann für d wegen der Spline- Definition d(j) (xi ) = 0, j = 0,... , m 1. Dies ergibt die Gestalt d(x) = bi (x xi ) m +, x 2 [xi , xi+1 ]. Wegen (x xi ) m + = 0 auf [x0 , xi ] ist die Behauptung auch für i + 1 richtig. 2. Die Linearität des Raumes Snm ([a, b]) ist o↵ensichtlich. Schließlich ergibt sich die Dimension aus Lemma 2.43 und der Analyse aus Teil (i) des Beweises. 2.3. SPLINE-INTERPOLATION 39 2.3.2 Interpolation in Spline-Räumen Bei Interpolation einer Funktion f durch Splines m ter Ordnung an den n+1 Stützstellen x0 ,..., xn , d.h. s(xj ) = f (xj ), j = 0,..., n, werden m 1 Freiheitsgrade des Raumes Snm ([a, b]) nicht genutzt. Eine Ausnahme bildet hier der Fall stückweise linearer Splines (m = 1, vgl. Beispiel 2.41). Die Menge Sn1 ([a, b]) ist unisolvent bezüglich der Stützstellen x0 ,... , xn. (Der Nachweis sei als Übungsaufgabe überlassen.) Für Splines mit m 2 kann man zusätzlich, zum Beispiel am Rand des Intervalls [a, b], Bedingungen stellen. Aus Symmetriegründen betrachten wir nur ungerade Zahlen m = 2l 1 mit l 2 N, l 2 für eine der Randbedingungen s(j) (a) = f (j) (a), s(j) (b) = f (j) (b), j = 1,... , l 1, (2.39a) bzw. s(l+j) (a) = s(l+j) (b) = 0, j = 0,... , l 2, (2.39b) oder für periodische Funktionen f mit Periode b a die Periodizitätsbedingungen2 s(j) (a) = s(j) (b), j = 1,... , m 1 = 2l 2, (2.39c) Im Fall (2.39a) spricht man auch von Hermite-Randbedingungen. Die sogenannten natürlichen Rand- bedingungen (2.39b) sind jedoch nicht von praktischer Bedeutung. Zunächst zeigen wir eine Extremaleigenschaft für Splines. Lemma 2.45 Sei f 2 C l [a, b] für l 2 N, l 2 und s 2 Snm ([a, b]) mit m = 2l 1 die interpolierende Spline- Funktion, d.h. s(xj ) = f (xj ), j = 0,..., n. Ferner gelte eine der Randbedingungen (2.39a), (2.39b) bzw. (2.39c). Dann gilt Z bh i2 Z bh i2 Z bh i2 f (l) (x) s(l) (x) dx = f (l) (x) dx s(l) (x) dx. a a a Beweis. Ausmultiplikation ergibt Z bh i2 Z bh i2 Z bh i2 f (l) (x) s(l) (x) dx = f (l) (x) dx s(l) (x) dx 2S a a a mit Z bh i S := f (l) (x) s(l) (x) s(l) (x)dx. a Auf S wenden wir (l 1)-fach die Regel der partiellen Integration an und berücksichtigen jeweils einer der Bedingungen (2.39a), (2.39b) bzw. (2.39c). Dann erhalten wir Z b⇥ ⇤ S = ( 1) l 1 f 0 (x) s0 (x) s(m) (x)dx, a 2 Die Bedingungen (2.39c) können auch für j = 0 verwendet werden, wenn die Stützstellen xn nicht als Interpolati- onsbedingung erfasst wird (x0 wird ja bereits erfasst). 40 KAPITEL 2. INTERPOLATION da f 2 C l [a, b] und für s 2 C m 1 [a, b] die Ableitung s(m) stückweise existiert. Erneute partielle Integration liefert mit den Interpolationsbedingungen sowie wegen s 2 ⇧m ([xj 1 , xj ]), dass n Z X xj ⇥ ⇤ S = ( 1)l 1 f 0 (x) s0 (x) s(m) (x)dx j=1 xj 1 n Z ! X x xj = ( 1)l 1 [f (x) s(x)] s(m) (x)|xjj 1 [f (x) s(x)] s(m+1) (x) dx = 0. j=1 xj 1 Daraus folgt die Behauptung. Lemma 2.45 ergibt die Ungleichung Z bh i2 Z bh i2 (l) s (x) dx  f (l) (x) dx. a a Dies erlaubt für den Fall kubischer Splines, d.h. m = 3, eine geometrische Interpretation. Für die Krümmung einer durch y = g(x) beschriebenen Kurve gilt g 00 (x) k(x) = , (1 + [g 0 (x)]2 )3/2 d.h. im Falle |g 0 (x)| ⌧ 1 ist näherungsweise k(x) ⇡ g 00 (x) bzw. Z b kkk22 := [k(x)]2 dx ⇡ kgk22. a Die abgeleitete Ungleichung bedeutet nun, dass der interpolierende kubische Spline unter allen Funk- tionen g 2 C 2 [a, b], die die gleichen Interpolationsforderungen erfüllen, näherungsweise die mittlere Krümmung kkk2 minimiert. Als für uns wesentlichste Folgerung aus Lemma 2.45 zeigen wir jetzt, dass die Bestimmung der in- terpolierenden Spline-Funktion s 2 Snm ([a, b]) mit m = 2l 1, l 1 an eine gegebene Funktion f unter einer der Randbedingungen (2.39a), (2.39b) oder (2.39c) wohldefiniert ist. Speziell muss man sich nicht explizit um die Stetigkeitsbedingungen in den Ableitungen s(j) für j = 1,... , 2l 2 der Spline-Funktion in den inneren Stützstellen x1 ,... , xn 1 kümmern. Genauer gilt Satz 2.46 Unter den Voraussetzungen von Lemma 2.45 existiert genau ein Spline s 2 Snm ([a, b]), der die Funktion f an den Stützstellen x0 ,..., xn interpoliert und eine der Randbedingungen (2.39a), (2.39b) oder (2.39c) erfüllt. Beweis. Unter Verwendung der Darstellung (2.38) durch Kardinalsplines ergibt sich aus den Inter- polationsforderungen und den Randbedingungen (2.39a), (2.39b) bzw. (2.39c) ein Gleichungssystem für die Koeffizienten ak und bk. Nach dem Beweis zum Satz 2.24 ist lediglich zu zeigen, dass der interpolierende Spline s zu f ⌘ 0 auf [a, b] identisch verschwindet. Rb⇥ ⇤2 Lemma 2.45 mit f ⌘ 0 liefert a s(l) (x) dx = 0, damit s(l) (x) = 0 bzw. s 2 ⇧l 1 [a, b]. Wegen der Randbedingungen für s folgt aber s(x) ⌘ 0. 2.3. SPLINE-INTERPOLATION 41 2.3.3 B-Splines Eine naheliegende Variante zur Berechnung der Koeffizienten ak , bk im Ansatz (2.38) wäre, das entstehende Gleichungssystem aus den Interpolationsforderungen s(xj ) = f (xj ), j = 0,... , n und einer der Randbedingungen (2.39a) oder (2.39c) zu lösen. Dieser Weg ist aber bei großer Zahl n 1 von Stützstellen nicht geeignet, da die entsprechende Koeffizientenmatrix sehr stark besetzt und schlecht konditioniert ist. Dies liegt vor allem daran, dass ein Teil der Kardinalsplines das gesamte Intervall [a, b] als Träger hat. Einen Ausweg stellt die Verwendung von Basisfunktionen in Snm ([a, b]) mit kleinem Träger dar, d.h. sie sind nur auf einem kleinem Teilgebiet des Intervalls [a, b] von Null verschieden. Auf Schönberg gehen die sogenannten B-Splines (basic spline curves) zurück, die wir vereinfachend für den äquidistanten Fall beschreiben. Definition 2.47 Ausgehend von ⇢ 1, |x|  12 , B0 (x) := 0, |x| > 12 definiert man B-Splines rekursiv durch Z x+ 1 2 Bm+1 (x) := Bm (t) dt, x 2 R, m = 0, 1,.... (2.40) 1 x 2 Durch vollständige Induktion nach m gewinnt man folgendes Resultat. Lemma 2.48 Die durch (2.40) definierten B-Splines Bm (·), m 2 N0 gehören zum Raum C m 1 (R), sind nichtne- gativ und verschwinden außerhalb des Intervalls [ m 2 1 m 1 2 , 2 + 2 ]. Sie sind stückweise polynomial vom Grad m auf den Intervallen [i, i + 1] für ungerade Zahlen m bzw. [i 12 , i + 12 ] für gerade Zahlen m für beliebige ganze Zahlen i. Beweis. (Übungsaufgabe !) Beispiel 2.49 (B-Spline der Ordnung m = 1, 2, 3) (i) Für m = 1 erhalten wir 8 < 1 |x|, |x|  1, B1 (x) = : 0, |x| 1. Diese wegen ihrer Form Hutfunktionen (hat functions) genannten Splines spielen bei der Näherungslösung von Di↵erential- und Integralgleichungen eine wichtige Rolle. Sie stellen das einfachste Beispiel für sogenannte finite Elemente dar. Man kann diese stückweise linearen, jedoch global stetigen Funktionen leicht auf den multiva- riaten Fall verallgemeinern, wenn man sie über simplizialen Gebieten im Rn definiert. Das sind im Fall n = 1 gerade Intervalle, für n = 2 Dreiecke und für n = 3 Tetraeder. Eine weitere Modifikation im univariaten Fall mit m = 1 ist, dass man global stetige, jedoch stückweise polynomiale Funktionen höherer Ordnung über den Intervallen [xj 1 , xj ] konstru- iert. Wir werden diese Idee bei der numerischen Integration in Kapitel 3.1 benutzen. 42 KAPITEL 2. INTERPOLATION (ii) Für m = 2 erhält man 8 1 2 > > 2 (|x| 2) (|x| + 12 )2 , |x|  12 , > > 1< 3 2 1 B2 (x) = (|x| 2) , 2  |x|  32 , 2> > > > : 3 0, |x| 2. (iii) Schließlich gewinnt man für m = 3 die Darstellung 8 > > (2 |x|)3 4(1 |x|)3 , |x|  1, > > 1< B3 (x) = (2 |x|)3 , 1  |x|  2, 6> > > > : 0, |x| 2. Abbildung 2.10: B-Splines der Ordnung m 2 {0, 1, 2, 3} Wir wollen jetzt bei äquidistanter Zerlegung eines Intervalles [a, b] zeigen, dass durch B-Splines eine Basis des Raumes Snm ([a, b]) erzeugt wird. Satz 2.50 Sei durch xk = a + hk mit k = 0,... , n und n 2 eine äquidistante Zerlegung des Intervalls [a, b] mit der Schrittweite h = n1 (b a) gegeben. Gelte ferner m = 2l 1 mit l 2 N. Dann erhält man eine Basis des Raumes Snm ([a, b]) durch die transformierten B-Splines ✓ ◆ x xk Bm,k (x) := Bm , k= l + 1,... , n + l 1, x 2 [a, b]. (2.41) h 2.3. SPLINE-INTERPOLATION 43 Beweis. (i) Wir beweisen zunächst durch Induktion nach m 2 N0 , dass die verschobenen B-Splines Bm (· k), k = 0,... , m ⇥1 ⇤ auf dem Intervall Im = 2 (m 1), 12 (m + 1) linear unabhängig sind. Für m = 0 ist die Behauptung nach Definition 2.47 o↵enbar erfüllt. Die Behauptung sei auch richtig für eine Zahl m 1 mit m 2 N, d.h. m 1 X  1 1 ↵k Bm 1 (x k) = 0, x 2 (m 2), m =) ↵k = 0, k = 0,... , m 1 2 2 k=0 Xm ✓ ◆  1 1 1 ↵k 1 Bm 1 x k+ = 0, x 2 (m 1), (m + 1) ) ↵k = 0, k = 0,..., m 1. 2 2 2 k=1 | {z } Im (2.42) Gelte nun m X k Bm (x k) = 0, x 2 Im. (2.43) k=0 Di↵erentiation in (2.43) ergibt unter Beachtung der Definition 2.47 von B-Splines X m  ✓ ◆ ✓ ◆ 1 1 k Bm 1 x k+ Bm 1 x k = 0, x 2 Im. 2 2 k=0 Es gilt nach Lemma 2.48 ✓ ✓ ◆◆ ✓ ◆ 1 m+1 m 1 supp Bm 1 ·+ = , , 2 2 2 ✓ ✓ ◆◆ ✓ ◆ 1 m + 1 3m + 1 und supp Bm 1 · m = ,. 2 2 2 Die Funktionen Bm 1 (· + 12 ) und Bm 1 (· m 1 2) verschwinden demnach auf Im. Daher gilt nach Indexverschiebung m X ✓ ◆ 1 ( k k 1 ) Bm 1 x k+ = 0, x 2 Im. 2 k=1 Damit folgt nach Induktionsvoraussetzung (2.42), dass k = k 1 für k = 1,... , m bzw. k = 0 = für k = 0,... , m. Somit können wir Formel (2.43) auch schreiben als m X Bm (x k) = 0, x 2 Im. k=0 Nach Integration über das Intervall Im ergibt sich Xm Z Xm Z m+1 k Z 1 (m+1) 2 2 Bm (x k) dx = Bm (t) dt = Bm (t) dt = 0, m 1 1 k=0 Im k=0 2 k 2 (m+1) | {z } >0 wo bei wir die Nichtnegativität der Splines Bm nach Lemma 2.48 ausgenutzt haben. Damit ist = 0. Daraus folgt die lineare Unabhängigkeit auch für die Zahl m und damit die Induktions- aussage. 44 KAPITEL 2. INTERPOLATION (ii) Man rechnet unmittelbar nach, dass die nach (2.41) definierten Splines ✓ ◆ ✓ ◆ x xk x a Bm,k (x) = Bm ⌘ Bm k , k = l + 1,... , n + l 1 h h unter Beachtung von Lemma 2.48 zum Raum Snm ([a, b]) gehören. Durch sinngemäße Anwen- dung der in Teil (i) des Beweises gezeigten Aussage folgt, dass diese m + n Funktionen linear unabhängig sind. Die Behauptung folgt dann aus Satz 2.44. Man kann die Spline-Interpolation mittels B-Spline-Basis wesentlich günstiger als mit Kardinal- Splines berechnen. Dies liegt wesentlich daran, dass die B-Splines nur einen kleinen Träger haben. Im Ergebnis erhält man eine im allgemeinen Fall große, jedoch schwachbesetzte Matrix. Wir demonstrieren dies am Beispiel kubischer Splines, d.h. für m = 3. Nach Beispiel 2.49 verschwindet B3 außerhalb des Intervalls ( 2, 2). Weiter ergibt sich 2 1 1 B3 (0) = , B3 (±1) = , B30 (0) = 0, B30 (±1) = ⌥. 3 6 2 Gesucht wird der kubische Spline n+1 X ✓ ◆ x xj s(x) = c j B3 , x 2 [a, b] (2.44) h j= 1 mit den Interpolationsbedingungen s(xj ) = yj := f (xj ) mit j = 0,... , n sowie den Randbedingun- gen s0 (a) = a1 , s0 (b) = b1. Die Randbedingungen ergeben c1 c 1 cn+1 cn 1 s0 (a) = = a1 , s0 (b) = = b1. 2h 2h Dann genügen die gesuchten Koeffizienten c 1 ,... , cn+1 dem folgenden linearen Gleichungssystem T AC = F, C = (c 1 ,..., cn+1 ) mit F = (ha1 , y0 , y1 ,... , yn , hb1 )T 2 Rn+3 und der schwach besetzten Matrix 0 1 3 0 3 B 1 4 1 C B C 1B B 1 4 1 C C A= B...... C 2 R(n+3)⇥(n+3). 6B... C B C @ 1 4 1 A 3 0 3 Die Matrix A erfüllt das schwache Zeilensummenkriterium und ist nicht reduzierbar (unzerlegbar). Dann existiert die inverse Matrix A 1. Ferner sind die grundlegenden Iterationsverfahren (siehe Ka- pitel [?]) zur Lösung des Gleichungssystems anwendbar. 2.3. SPLINE-INTERPOLATION 45 2.3.4 Fehlerabschätzungen für Splines Für den Fall stückweise linearer Splines aus Sn1 ([a, b]) hatten wir im Beispiel 2.41 bereits Fehler- abschätzungen hergeleitet. Entsprechende Untersuchungen für Spline-Räume Snm ([a, b]) mit m > 1 sind wesentlich aufwendiger. Wir zeigen hier exemplarisch für den Fall kubischer Splines Fehlerabschätzungen, die geeignet verall- gemeinert werden können. Zunächst betrachten wir Funktionen f 2 C 2 [a,