Interpolation I

Computer-Algebra und Visualisierungen von Interpolationsformeln

Prof. Dr. Bernhard Zgraggen, HSR

Einführung in dieses Notebook

Dieses Dokument ist ein "selbsterklärendes" Mathematica-Notebook, das interaktiv bearbeitet werden kann ("Hands-On"). Input-Zellen können mittels der Tastenkombination ⇧⌤ (in Worten: Shift-Enter) ausgeführt werden. Falls Textzellen geschlossen sind, so können sie durch Klick auf das Icon ▼ oder auf die blaue Zellklammer rechterhand geöffnet werden.

Wenn Sie in Mathematica auf Befehle oder Begriffe stossen, die Sie nicht verstehen, so steht Ihnen die mächtige Hilfe/Doc Center zur Verfügung; diese ist interaktiv.

Allfällige Aufgaben/Übungen befinden sich am Ende des Dokuments.

Symbolische Daten, Newton-Basispolynome, Dividierte Differenzen → {InterpolatingPolynomial}

Interpolationspolynom ohne wiederholte Argumente (ohne Ableitungen = Collocation)

Definition symbolischer Datenpaare in einer Liste variabler Länge.

datasymb[n_]:=Table[{x[i],y[i]},{i,0,n}]

datasymb[4]//TableForm

Definition der Newton-Polynome.

newtonpoly[k_]:=Product[#-x[i],{i,0,k-1}]&

Table[newtonpoly[k][x],{k,0,4}]//TableForm

Auswertung der Newton-Polynome an den Stellen x[0], x[1], ... und Gleichungssystem der Interpolation:

Manipulate[StringForm["`1` == `2`",TableForm@Table[(∑_(k=0)^n)a[k](newtonpoly[k][x[i]]),{i,0,n}],TableForm@Table[y[i],{i,0,n}]],
{{n,0,"Anzahl Argumente - 1"},0,10,1,Appearance->"Labeled"}
]

Die Berechnung eines symbolischen Interpolationspolynoms zeigt den Aufbau mit dividierten Differenzen (hier in der HornerForm ; diese ist effizient und numerisch eher stabil für Auswertungen).

InterpolatingPolynomial[datasymb[4],x]

Das Ausmultiplizieren (Expand) unter Beibehaltung der Differenzen (x - x[i]) führt zur Lagrange-Darstellung des Interpolationspolynoms.

Expand[%,y[_]];

Eine Art interaktives "Instant GUI" kann mit Manipulate erzeugt werden. Dieser Befehl ist sehr vielseitig und mächtig.

Manipulate[InterpolatingPolynomial[datasymb[n],x],{{n,0,"Anzahl Argumente - 1"},0,10,1,Appearance->"Labeled"}]

Die Rekursionsbeziehung von Neville-Aitken (Collocation)

Der folgende Code zeigt die Korrektheit der Beziehung von Neville-Aitken.

n=4;
(InterpolatingPolynomial[datasymb[n],x]
-((1/((x[n]-x[0])))((x-x[0])InterpolatingPolynomial[Rest@datasymb[n],x]-(x-x[n])InterpolatingPolynomial[Most@datasymb[n],x])))//Simplify

Das Schema der dividierten Differenzen (erweitert um gleichverteilte Argumente)

Mit der folgenden Code-Sequenz wird das Schema der dividierten Differenzen rekursiv definiert. Die Definition ist rekursiv und bedingt (Conditional-Operator /; ). Die x-Argumente werden als Liste übergeben, bspw. {x[0], x[1], x[2]}. Wenn man die Übergaben mit Liste vermeiden möchte könnte man den Sequence Operator verwenden.

Clear[x,y]
y[{x[i_Integer]}]:=y[i]

y[args_List]:=((y[Rest@args]-y[Most@args])/(Last[args]-First[args]))/;(Length[args]≥2 )

y[{x[0],x[1],x[2],x[3]}]

n=3;
Table[Table[y[Table[x[j],{j,i,k}]],{i,k,0,-1}],{k,0,n}]//TableForm

Eine Art interaktives "Instant GUI" kann mit Manipulate erzeugt werden. Dieser Befehl ist sehr vielseitig und mächtig.

Manipulate[Table[Table[y[Table[x[j],{j,i,k}]],{i,k,0,-1}],{k,0,nn}]//TableForm,{{nn,0,"Anzahl Argumente - 1"},0,10,1,Appearance->"Labeled"}]

Im folgenden die Verwendung numerischer Daten mittels ReplaceAll /. . Zuweisungen an x[0] etc. und y[0] etc. funktionieren hier nicht (Konflikt mit der Definition von y oben.)

datanum={x[0]->0,x[1]->2,x[2]->5,x[3]->-2,x[4]->-1,y[0]->-1, y[1]->1, y[2]->2,y[3]->0,y[4]->1}

n=4;
(Table[Table[y[Table[x[j],{j,i,k}]],{i,k,0,-1}],{k,0,n}]/.datanum)//TableForm

Numerische Daten lassen sich wie folgt mit Manipulate verwenden.

datanum={x[0]->0,x[1]->2,x[2]->5,x[3]->-2,x[4]->-1,y[0]->-1, y[1]->1, y[2]->2,y[3]->0,y[4]->1}

n=4;
Manipulate[(Table[Table[y[Table[x[j],{j,i,k}]],{i,k,0,-1}],{k,0,nn}]/.datanum)//TableForm,{{nn,0,"Anzahl Argumente - 1"},0,n,1,Appearance->"Labeled"}]

Bei gleichmässigen Argumenten ergibt sich der Spezialfall mit gewöhnlichen Differenzen:

n=4;
datanum=Table[x[i]->x[0]+h i,{i,0,n}]

Manipulate[(Table[Table[y[Table[x[j],{j,i,k}]],{i,k,0,-1}],{k,0,nn}]/.datanum)//TableForm,{{nn,0,"Anzahl Argumente - 1"},0,n,1,Appearance->"Labeled"}]

Interpolationspolynom mit wiederholten Argumente (inkl. Ableitungen = Osculation)

Definition symbolischer Datenpaare in einer Liste variabler Länge inkl. Ableitungen höherer Ordnung. Die Ordnungen der Ableitungen sind in der Liste orders gespeichert. Die Länge dieser Liste ist gleich n+1 und ihre Einträge sind natürliche Zahlen.

datasymb[n_,orders_]:=Table[Flatten[{x[i],y[i],Table[(Subscript[y, j])[i],{j,1,orders[[i+1]]}]}],{i,0,n}]

Hier werden die Ableitungen der Ordnungen 1, 1, 1, 3, 2 für die sukzessiven x-Argumente vorgegeben.

datasymb[4,{1,1,1,3,2}]//TableForm

Die Berechnung eines symbolischen Interpolationspolynoms zeigt den Aufbau mit dividierten Differenzen (hier in der HornerForm ; diese ist effizient und numerisch eher stabil für Auswertungen).

datasymb[2,{0,1,0}]//TableForm

InterpolatingPolynomial[datasymb[2,{0,1,0}],x]

Das Ausmultiplizieren (Expand) unter Beibehaltung der Differenzen (x - x[i]) führt zur Lagrange-Darstellung des Interpolationspolynoms.

Expand[%,y[_]]

Eine Art interaktives "Instant GUI" kann mit Manipulate erzeugt werden. Dieser Befehl ist sehr vielseitig und mächtig. Die Ableitungsordnung ist für alle Argumente gleich.

Manipulate[InterpolatingPolynomial[datasymb[n,Table[order,{n+1}]],x],
{{n,0,"Anzahl Argumente - 1"},0,10,1,Appearance->"Labeled"},{{order,0,"Ableitungsordnung"},0,10,1,Appearance->"Labeled"}]

Interpolationsfehler

Der Interpolationsfehler ohne wiederholte Argumente (Collocation)

Die grundlegende Formel für den Interpolationsfehler lässt sich für geeignete theoretische Funktionen wie folgt ausdrücken:

   y(x)-p(x)=(((y^((n+1)))(ξ))/((n+1)!))(Subscript[ℼ, n+1])(x)=(((y^((n+1)))(ξ))/((n+1)!))(x-(Subscript[x, 0]))(x-(Subscript[x, 1]))⋯(x-(Subscript[x, n-1]))(x-(Subscript[x, n]))         ξ, x ∈ (min{(Subscript[x, 0]),} max{(Subscript[x, i])})

Hierbei bezeichnen n+1 die Anzahl Argumente ((Subscript[x, i])) bzw. n den Grad des Interpolationspolynoms. Das Verhalten dieses Fehlers in Abhängigkeit von n und x ist abhängig vom Verhalten der Ableitung (((y^((n+1)))(ξ))/((n+1)!)) und des Produkts der Differenzen (Subscript[ℼ, n+1])(x). Ein einfaches Beispiel für Fehlerdivergenz ist das Runge-Phänomen.

Der Interpolationsfehler bei wiederholten Argumenten (Vorgabe von Ableitungswerten = Osculation)

In Analogie zur grundlegenden Formel für den Interpolationsfehler bei der Collocation lässt sich für geeignete theoretische Funktionen der Fehler bei Osculation wie folgt ausdrücken. Hierbei bezeichnen a die Gesamtanzahl Bedingungen an die Daten sowie (Subscript[a, i]) die Anzahl Bedingungen an das Argument (Subscript[x, i]).

   y(x)-p(x)=(((y^((a)))(ξ))/((a)!))(Subscript[ℼ, n+1])(x)=(((y^((n+1)))(ξ))/((n+1)!))(((x-(Subscript[x, 0])))^(Subscript[a, 0]))(((x-(Subscript[x, 1])))^(Subscript[a, 1]))⋯(((x-(Subscript[x, n-1])))^(Subscript[a, n-1]))(((x-(Subscript[x, n])))^(Subscript[a, n]))        ξ, x ∈ (min{(Subscript[x, 0]),} max{(Subscript[x, i])})

Das Runge-Phänomen

In der folgenden Zelle werden die Funktionsformel y (zum Beispiel initial die Runge-Funktion), die Argumente und ein Intervall für den Plotbereich (x-Achse) definiert. Danach kann generisch ein Plot erzeugt werden. Darunter befindet sich zusätzlich ein generisches Manipulate ("Instant GUI").

y[x_]:=1./(1+25x^2);
n=25; (* Anzahl Argumente - 1 = Grad des Polynoms *)
{a,b}={-1.,1.};

datapts=Table[{x,y[x]},{x,a,b,(b-a)/(n)}]; (* Datenpunkte als Paare {x, y} *)
GraphicsRow[
{Plot[Evaluate@InterpolatingPolynomial[datapts,x],{x,a,b},PlotRange->{{a,b},Automatic},Epilog->{Red,Map[Point,datapts]}, PlotLabel->"Interpolation"],
Plot[Evaluate@Product[x-datapts[[i+1]][[1]],{i,0,n}],{x,a,b},PlotRange->{{a,b},Automatic},PlotStyle->Darker[Red],PlotLabel->"Newtonprodukt (Subscript[ℼ, n + 1])(x)"],
Plot[Evaluate@(D[y[x],{x,n+1}]/(n+1)!),{x,a,b},PlotRange->{{a,b},Automatic},PlotStyle->Darker[Blue],PlotLabel->"Ableitungsterm"]
}
]

y[x_]:=1./(1+25x^2);
n=25; (* Anzahl Argumente - 1 = Grad des Polynoms *)
{a,b}={-1,1};

demoRunge1=Manipulate[Module[{datapts},
If[!rand&&!cheb,
datapts=Table[{x,y[x]},{x,a,b,(b-a)/(n)}] (* Datenpunkte als Paare {x, y} *),
If[rand,datapts={#,y[#]}& /@RandomReal[{a,b},{n+1}]];
If[cheb,datapts={#,y[#]}& /@Table[a+((b-a)/2)(Cos[(((2i+1)Pi)/(2(n+1)))]+1),{i,0,n}]]
]; 
GraphicsRow[
{Plot[{Evaluate@InterpolatingPolynomial[datapts,x],Evaluate@(y[x])},{x,a,b},PlotRange->{{a,b},Automatic},PlotStyle->{Automatic,{Red,Dashed}},Epilog->{Red,Map[Point,datapts]}, PlotLabel->"Interpolation"],
Plot[Evaluate@Product[x-datapts[[i+1]][[1]],{i,0,n}],{x,a,b},PlotRange->{{a,b},Automatic},PlotStyle->Darker[Red],PlotLabel->"Newtonprodukt (Subscript[ℼ, n + 1])(x)"],
Plot[Evaluate@(D[y[x],{x,n+1}]/((n+1)!)),{x,a,b},PlotRange->{{a, b},Automatic},PlotStyle->Darker[Blue],PlotLabel->"Ableitungsterm"]
}
]
] (* end of Module *),
{{n,10,"n"},1,100,1,Appearance->"Labeled"},
{{rand,False,"random args"},{True,False}},
{{cheb,False,"Chebyshev args"},{True,False}}] (* end of Manipulate *)

Mehrdimensionale Interpolationen

Erklärung mit Visualisierungen (Collocation)

Die Interpolation mehrdimensionaler Daten kann stufenweise auf tieferdimensionale Interpolationen zurückgeführt werden. Besonders häufig sind die Fälle 2D (zum Beispiel Bildverarbeitung) und 3D (zum Beispiel 3D-Grafikprogrammierung) und damit verbunden Begriffe wie Bi-Linear/Bi-Kubisch bzw. Tri-Linear/Tri-Kubisch etc.

Das grundsätzliche Vorgehen wird am 2D-Fall illustriert. Es soll für eine Menge M von Argumentpaaren {x,y} ein Interpolations-Polynom p(x,y) in 2 Variablen x und y gefunden werden, das bei den vorgegebenen Argumentpaaren aus M vorgegebene z-Werte annimmt. Eine solche Kollokations-Aufgabe lässt sich nur dann allgemein lösen, wenn die Daten genügend Freiheitsgrade aufweisen (ein typisches Gegenbeispiel basiert auf drei auf einer Geraden liegenden Argumentpaaren mit mind. zwei verschiedenen z-Werten).

Ein häufiger Fall für die Menge M der Punktepaare ist ein Gitter (grid): {x[0], x[1], ..., x[m]} × {y[0], y[1], ...,  y[n]}.

Die folgende Sequenz definiert ein Gitter von ganzzahligen Argumenten sowie zugehörige z-Werte.

{m,n}={3,5}; (* Gitterlängen *)
{xs,ys}={Range[0,m],Range[0,n]}; (* x-Werte, y-Werte *)
zs=RandomReal[{-1.,1.},{m+1,n+1}]; 
(* z-Werte *)

Daraus werden 3D-Punkte (Koordinatentripel) generiert und anschliessend illustriert.

datapts3d=Flatten/@Transpose[{Tuples[{xs,ys}],Flatten@zs}];

lp3d=ListPointPlot3D[datapts3d,Filling->Axis,FillingStyle->Dashed,Axes->True,AxesLabel->{"x","y","z"},PlotStyle->{ PointSize[Medium],Blue},AspectRatio->Automatic]

Interpolationspolynome in x-Richtung (Variable x) für die vorgegebenen y-Werte (Liste von Polynomen in x, die abhängen von den z-Werten bei den y-Gitterargumenten).

xinterpoly=Map[InterpolatingPolynomial[Transpose[{xs,#}],x]&,Transpose@zs];

Diese Polynomschar zusammen mit den 3D-Datenpunkten werden illustriert:

Show@{ParametricPlot3D[Table[{x,ys[[i]],xinterpoly[[i]]},{i,1,n+1}],{x,0,m},BoxRatios->{1,1,0.5},Axes->True,AxesLabel->{"x","y","z"}],lp3d}

In diese Polynomfunktionen lassen sich beliebige x-Werte einsetzen, entsprechende Datenpunkte berechnen und bzgl. der Variable y interpolieren. Der resultierende Ausdruck ist ein Polynom in den freien Variablen x und y.

yxinterpoly=InterpolatingPolynomial[Transpose[{ys,xinterpoly}],y];
%//Simplify

Der folgende Plot zeigt dieses Polynom als Flächenfunktion über der Koordinatenebene sowie ein Interpolationspolynom in der Variable y für einen beliebigen x-Wert.

x0=RandomReal[{0,m}];
Show[{
ParametricPlot3D[Table[{x,ys[[i]],xinterpoly[[i]]},{i,1,n+1}],{x,0,m},BoxRatios->{1,1,0.5},Axes->True,AxesLabel->{"x","y","z"}],ParametricPlot3D[Evaluate@{x0,y,yxinterpoly/.{x->x0}},{y,0,n},PlotStyle->{Darker[Red],Thick}],
Graphics3D[{PointSize[0.015],Red,Point[Table[{x0,ys[[i]],xinterpoly[[i]]/.{x->x0}},{i,1,n+1}]]}],
Plot3D[yxinterpoly,{x,0,m},{y,0,n},PlotStyle->Opacity[0.5],Mesh->Automatic,MeshStyle->Dotted,ColorFunction->Hue],lp3d},PlotRange->Automatic,BoxRatios->{1,1,0.5},Axes->True,AxesLabel->{"x","y","z"}]

Der folgende ContourPlot zeigt dieses Polynom in Konturendarstellung über der Koordinatenebene. Die Farbfunktion ist durch die Transformation der z-Koordinate in das Intervall [0,1] und den entsprechenden Hue-Wert bestimmt.

ContourPlot[yxinterpoly,{x,0,m},{y,0,n},ColorFunction->Hue,Contours->128]

Integrierte dynamische Visualisierung (Collocation)

{m,n}={3,3}; (* Gitterlängen *)
{xs,ys}={Range[0,m],Range[0,n]}; (* x-Werte, y-Werte *)
zs=RandomReal[{-1.,1.},{m+1,n+1}]; 
(* z-Werte *)

demo1:=Manipulate[Module[{datapts3d,lp3d,xinterpoly,yxinterpoly,x,y},
datapts3d=Flatten/@Transpose[{Tuples[{xs,ys}],Flatten@zs}];
lp3d=ListPointPlot3D[datapts3d,Filling->Axis,FillingStyle->Dashed,Axes->True,AxesLabel->{"x","y","z"},PlotStyle->{PointSize[ps],Darker[Blue]}];
xinterpoly=Map[InterpolatingPolynomial[Transpose[{xs,#}],x]&,Transpose@zs];
yxinterpoly=InterpolatingPolynomial[Transpose[{ys,xinterpoly}],y];
Show[{
ParametricPlot3D[Table[{x,ys[[i]],xinterpoly[[i]]},{i,1,n+1}],{x,0,m},BoxRatios->{1,1,0.5},Axes->True,AxesLabel->{"x","y","z"}],
Graphics3D[{PointSize[ps],Darker[Red],Point[Table[{x0,ys[[i]],xinterpoly[[i]]/.{x->x0}},{i,1,n+1}]]}],
ParametricPlot3D[Evaluate@{x0,y,yxinterpoly/.{x->x0}},{y,0,n},PlotStyle->{Darker[Red],Thick}],
ParametricPlot3D[Evaluate@{x,y0,yxinterpoly/.{y->y0}},{x,0,m},PlotStyle->{Red,Thick}],
Plot3D[yxinterpoly,{x,0,m},{y,0,n},PlotStyle->Opacity[0.5],Mesh->True,MeshStyle->{{Dotted,Gray}},ColorFunction->Hue],
lp3d},
PlotRange->Automatic,BoxRatios->{1,1,0.5},Axes->True,AxesLabel->{"x","y","z"}]
]
,
{{x0,Min[xs],"x0"},Min[xs],Max[xs],Appearance->"Labeled"},
{{y0,Min[ys],"y0"},Min[ys],Max[ys],Appearance->"Labeled"},
{{ps,0.01,"Rel. Punktgrösse für Schnitt"},0.,0.1,Appearance->"Labeled"}]

demo1

Rein symbolische Daten (zur Formel-Erklärung)

Die folgende Sequenz definiert ein Gitter von symbolischen Argumenten sowie zugehörige z-Werte.

{m,n}={2,3}; (* Gitterlängen *)
{xs,ys}={Table[x[i],{i,0,m}],Table[y[i],{i,0,n}]}; (* x-Werte, y-Werte *)
zs=Table[z[i,j],{i,0,m},{j,0,n}];

Daraus werden 3D-Punkte (Koordinatentripel) generiert:

datapts3d=Flatten/@Transpose[{Tuples[{xs,ys}],Flatten@zs}]

Interpolationspolynome in x-Richtung (Variable x) für die vorgegebenen y-Werte  (Liste von Polynomen in x, die abhängen von den z-Werten bei den y-Gitterargumenten).

xinterpoly=Map[InterpolatingPolynomial[Transpose[{xs,#}],x]&,Transpose@zs]

In diese Polynomfunktionen lassen sich beliebige x-Werte einsetzen, entsprechende Datenpunkte berechnen und bzgl. der Variable y interpolieren. Der resultierende Ausdruck ist ein Polynom in den freien Variablen x und y.

yxinterpoly=InterpolatingPolynomial[Transpose[{ys,xinterpoly}],y];

Expand[%,z[_,_]];

Bilinear-Interpolation auf Einheitsquadrat mit ContourPlot

{m,n}={1,1}; (* Gitterlängen *)
{xs,ys}={Range[0,m],Range[0,n]}; (* x-Werte, y-Werte *)
zs=RandomReal[{-1.,1.},{m+1,n+1}]; 
(* z-Werte *)

(* Script-Example *)
zs={{91,162},{210,95}}/255.
(* z-Werte deterministisch*)

datapts=Transpose[{Tuples[{xs,ys}],Flatten@zs}]

ip2d=InterpolatingPolynomial[datapts,{x,y}]

%//Expand

ContourPlot[ip2d,{x,0,m},{y,0,n},ColorFunction->Function[{z},RGBColor[z,0,0]],ColorFunctionScaling->True,Contours->32]

demo1

Trilinear-Interpolation auf Einheitswürfel mit ContourPlot3D

Image3D[{,,}]

{n1,n2,n3}={1,1,1}; (* Gitterlängen *)
{x1s,x2s,x3s}={Range[0,n1],Range[0,n2],Range[0,n3]}; (* x1-Werte ... *)
zs=RandomReal[{-1.,1.},{n1+1,n2+1,n3+1}];
(* z-Werte *)

(* z-Werte deterministisch*)

datapts=Transpose[{Tuples[{x1s,x2s,x3s}],Flatten@zs}]

ip3d=Interpolation[datapts,{x1,x2,x3}]

ContourPlot3D[ip3d,{x1,0,n1},{x2,0,n2},{x3,0,n3},ColorFunction->Hue,Contours->2*8]

*Aufgabe(n)

(1)
Mit den Tools für die Schemen der dividierten Differenzen und den Interpolationspolynomen können Sie die Übungen zu den Handrechnungen kontrollieren und dynamisch visualisieren.

(2)
Untersuchen Sie mit den Tools in diesem Notebook, ob es ein Runge-Phänomen gibt für die folgenden Modellfunktionen:

a) y = Log[1+x]    (-0.999999 ≤ x ≤ 0.999999)            b) y = Sin[x]/x    (-2.ℼ ≤ x ≤ 2.ℼ und für x = 0 ist der Wert y = 1.)      

c) y = Exp(-1/x^2)   (-1≤ x ≤ 1 und für x = 0 ist der Wert y = 0)

Untersuchen Sie jeweils zufällig verteilte Argumente und Chebyshev-Argumente.

Hinweise zu bc): Die Definition der Funktion y in b) erfolgt wie folgt: y[x_]:=Piecewise[{{1. Sin[x]/x, x ≠ 0}, {1., x == 0}}];

(3)
Führen Sie eine Bilinear-Interpolation von Hand durch für das Einheitsquadrat-Gitter und die z-Werte: 

{xs, ys} = {{0, 1}, {0, 1}}
zs = {{0, 1}, {1, 0.5}}

Visualisieren Sie das Ergebnis mittels ContourPlot.

Lösung: 

{xs, ys} = {{0, 1}, {0, 1}}
zs = {{0, 1}, {1, 0.5}}
datapts = Transpose[{Tuples[{xs, ys}], Flatten@zs}]
InterpolatingPolynomial[datapts, {x, y}]
ContourPlot[%,{x,0,1},{y,0,1}, ColorFunction->Hue,Contours->128]