Zusammenstellung ausgewählter Interpolations- und Approximationsverfahren
Lars Ulke-Winter
BTU Cottbus-Senftenberg
Fachgebiet Polymerbasierter Leichtbau
lars.ulke-winter@b-tu.de

Die beschriebenen Methoden und Built-In Subroutinen werden o. B. d. A. an Beispielen eindimensionaler Funktionen und generierten Datensätze dargestellt.

Clear["Global`*"]

Interpolationen

Ansatz nach Kriging

siehe auch: Veröffentlichung und Buch: Strukturoptimierung Abschnitt 6.2.4

(OverTilde[y])(x)=μ+rᵀ (R^(-1))((Subscript[y, i])-1μ)

Beispielfunktion mit Stützwerten

f[x_]:=Exp[x Cos[3π x]]
interval={x,0,1};

numberOfPkts=6;
xi=Table[(k/numberOfPkts),{k,0,numberOfPkts}]//N;
yi=f/@xi;

Abstandsfunktion

Θ=1; p=2;
d[xi_,xj_]:=Θ(((Abs[xi-xj]))^p)

Korrelationsmatrix R

rij=Table[Exp[-d[xi[[i]],xi[[j]]]],{i,numberOfPkts+1},{j,numberOfPkts+1}];
rij//MatrixForm

Korrelationsvektor r

ri=(Exp[-d[x,#]]&)/@xi

Mittelwert μ

μ=(ConstantArray[1,numberOfPkts+1]).Inverse[rij].yi/(ConstantArray[1,numberOfPkts+1]).Inverse[rij].(ConstantArray[1,numberOfPkts+1])

Interpolationsfunktion

fnKriging[x_]:=μ+ri.Inverse[rij].(yi-μ ConstantArray[1,numberOfPkts+1])

fnKriging[x]//Simplify

Plot[{f[x],fnKriging[x]},interval,PlotStyle->{Blue,Green},PlotLegends->{"Funktion","Interp. Kriging"},Epilog->Point[{xi,yi}\[Transpose]]]

Ansatz mit radialen Basisfunktionen
(RBF Radial Basis Functions)

ClearAll[f,interval,numberOfPkts,xi,yi,ϕ,m,c,fnRBF]

siehe auch Buch: Strukturoptimierung Abschnitt 6.2.4

Generiere Stützwerte aus Beispielfunktion

f[x_]:=Exp[x Cos[3π x]]
interval={x,0,1};

numberOfPkts=6;
xi=Table[(k/numberOfPkts),{k,0,numberOfPkts}]//N;
yi=f/@xi;

Basisfunktion deren Linearkombination die Funktion approximiert

ϕ[r_]:=Exp[-(((3 r))^2)]

(*ϕ[r_]:=r
ϕ[r_]:=(1/(Sqrt[1+((r^2)/(2^2))]))
ϕ[r_]:=(Sqrt[1+((r^2)/(2^2))])
ϕ[r_]:=(r^2)Log[r]*)

Koeffizientenmatrix erstellen

m=Table[ϕ[(((xi[[j]]-xi[[i]]))^2)],{i,numberOfPkts+1},{j,numberOfPkts+1}];

Zeile und Spalte mit 1 hinzufügen der Koeffizientmatrix hinzufügen.
Das n×n-Element der Koeffizientenmatrix ist 0:

AppendTo[m,ConstantArray[1,numberOfPkts+1]];
m=ResourceFunction["AppendColumn"][m,ConstantArray[1,numberOfPkts+2]];
m[[-1,-1]]=0;
m//MatrixForm

Ermittlung der Wichtungsfaktoren

c=LinearSolve[m,AppendTo[yi,0]]

Näherungsfunktion als Linearkombination und konstantem Glied

fnRBF[x_]:=(Table[ϕ[(((x-xi[[i]]))^2)],{i,numberOfPkts+1}].c[[1;;-2]])+c[[-1]]
fnRBF[x]//Simplify

Plot[{f[x],fnRBF[x]},interval,PlotStyle->{Blue,Green},PlotLegends->{"Funktion","Interp. RBF"},Epilog->Point[{xi,yi[[1;;-2]]}\[Transpose]]]

Interpolationen mit Wolfram Language Built-In Subroutinen

ClearAll[f,interval,numberOfPkts,xi,yi,fnPoly]

Generiere Stützwerte aus Beispielfunktion

f[x_]:=Exp[x Cos[3π x]]
interval={x,0,1};

numberOfPkts=6;
xi=Table[(k/numberOfPkts),{k,0,numberOfPkts}]//N;
yi=f/@xi;

Polynominterpolation

fnPoly[x_]:=InterpolatingPolynomial[{xi,yi}\[Transpose],x]//Expand
fnPoly[x]

Plot[{f[x],fnPoly[x]},interval,PlotStyle->{Blue,Green},PlotLegends->{"Funktion","Interp. Polynom"},Epilog->Point[{xi,yi}\[Transpose]]]

Das den n Stützstellen zugehörige Polynom vom Grade (n-1) ist eindeutig. 
Die Lösung eines zugehörigen Gleichungssystems kann vermieden werden, wenn sog. Lagrange-Polynome (Subscript[L, i])(x) mit der Eigenschaft (Subscript[L, i])((Subscript[x, k]))=({{{, ({{1,, für i=k}, {0,, für i≠k}})}})  aufgestellt werden; siehe auch: Lagrangesche Interpolationsformel und Buch Abschnitt A .1

Abschnittsweise Interpolation:
Vergleich Hermite- und Spline-Ansätze verschiedener Ordnung

(Plot[{f[x],Interpolation[{xi,yi}\[Transpose], Method->#[[1]],InterpolationOrder->#[[2]]][x]},interval,PlotStyle->{Blue,Green},PlotLabel->#[[3]], PlotLegends->{"Funktion","Interp."},Epilog->Point[{xi,yi}\[Transpose]]]&/@{{"Hermite",1,"Hermite 1st Order"},{"Hermite",2,"Hermite 2nd Order"},
{"Hermite",3,"Hermite 3rd Order"},
{"Spline",1,"Spline 1st Order"},
{"Spline",2,"Spline 2nd Order"},
{"Spline",3,"Spline 3rd Order"},
{"Spline",4,"Spline 4th Order"}
})//Column

Approximationen

Mit orthogonalen Ansatzfunktionen

ClearAll[f,interval,norm,n,ansatz,basePoly,baseTrig,baseRand,x,cPoly,cTrig,cRand,fPoly,fTrig,fRand]

Zu approximierende Beispielfunktion auf einem Intervall

f[x_]:=Exp[x Cos[3π x]]
interval={x,0,1};

Definition einer Norm

Norm auf einem Intervall (0,1):

norm[f_,g_]:=NIntegrate [f g,interval]

Orthogonalisieren linear unabhängiger Funktionen

Polynome:

n=6;
ansatz=Table[(x^i),{i,Range[0,n]}]

Orthogonalisierung (Orthonormalisierung) bspw. durch Gram Schmidt  Verfahren

basePoly=Orthogonalize[ansatz,Integrate[#1 #2,interval]&,Method->"GramSchmidt"];

Test der Orthonormalität (analog (Subscript[e, i])·(Subscript[e, j])=(Subscript[δ, ij])):

Integrate[basePoly[[1]] basePoly[[1]],interval]
Integrate[basePoly[[1]] basePoly[[2]],interval]

Trigonometrischen Funktionen (Fourier Reihe):

n=6;
ansatz=Prepend[Table[Sin[i x],{i,Range[1,n]}],1]

baseTrig=Orthogonalize[ansatz,NIntegrate[#1 #2,interval]&,Method->"GramSchmidt"]//Quiet;

Test der Ortnormalität (analog (Subscript[e, i])·(Subscript[e, j])=(Subscript[δ, ij])):

Integrate[baseTrig[[1]] baseTrig[[1]],interval]//Chop
Integrate[baseTrig[[1]] baseTrig[[6]],interval]//Chop

Menge wahllos zusammengestellter Ansatzfunktionen:

ansatz={1,x,Exp[x],Sin[x],Sin[2x],Log[x],Sqrt[x]}

baseRand=Orthogonalize[ansatz,NIntegrate[#1 #2,interval]&,Method->"GramSchmidt"]//Quiet;

NIntegrate[baseRand[[1]] baseRand[[1]],interval]//Quiet//Chop
NIntegrate[baseRand[[1]] baseRand[[2]],interval]//Quiet//Chop

Bestimmen der zugehörigen Wichtungsfaktoren

- Analog zum Skalarprodukt (Norm) der linearen Algebra: (Subscript[c, i]) (Subscript[e, i])·(Subscript[e, j])=(Subscript[c, i])(Subscript[δ, ij])=(Subscript[c, j])
- Siehe auch Galerkin Methode: Buch Abschnitt 3.6

cPoly=Table[norm[basePoly[[i]],f[x]],{i,Length[basePoly]}]//N

cTrig=Table[norm[baseTrig[[i]],f[x]],{i,Length[baseTrig]}]//Quiet//N

cRand=Table[norm[baseRand[[i]],f[x]],{i,Length[baseRand]}]//Quiet//N

Funktionsapproximation als Linearkombination orthonormaler Basisfunktionen

fPoly[x_]:=(∑_(i=1)^(Length[basePoly]))cPoly[[i]]basePoly[[i]]
fPoly[x]//Simplify

fTrig[x_]:=(∑_(i=1)^(Length[baseTrig]))cTrig[[i]]baseTrig[[i]]
fTrig[x]//Simplify

fRand[x_]:=(∑_(i=1)^(Length[baseRand]))cRand[[i]]baseRand[[i]]
fRand[x]//Simplify

Plot[{f[x],fPoly[x],fTrig[x],fRand[x]},interval,PlotStyle->{Blue,Green,Orange,Red},PlotLegends->{"Funktion","Approx. Poly.","Approx. Trig.","Approx. Rand."}]

Legendre Polynome:

ClearAll[interval,x,baseLegendre,cLegendre,fLegendre]

Legendre Polynome sind bereits orthogonal im Intervall (-1,1), jedoch i. Allg. nicht normiert

interval={x,-1,1};

baseLegendre=LegendreP[Range[0,12],x]

Test der Orthonormalität (analog (Subscript[e, i])·(Subscript[e, j])=(Subscript[δ, ij])):

Integrate[baseLegendre[[1]] baseLegendre[[1]],interval]//Chop
Integrate[baseLegendre[[1]] baseLegendre[[2]],interval]//Chop

Normalisierung der Basisfunktionen und Berechnung der Wichtungsfaktoren s.o.:

cLegendre=Table[((NIntegrate[baseLegendre[[i]]f[x],interval])/(NIntegrate[baseLegendre[[i]]baseLegendre[[i]],interval])) ,{i,Length[baseLegendre]}]//Quiet//N

Funktionsapproximation als Linearkombination orthonormaler Legendre Polynomen:

fLegendre[x_]:=(∑_(i=1)^(Length[baseLegendre]))cLegendre[[i]]baseLegendre[[i]]

fLegendre[x]//Simplify

Plot[{f[x],fLegendre[x]},interval,PlotStyle->{Blue,Green},PlotLegends->{"Funktion","Approx. Legendre"}]

Mit künstlichen Neuronalen Netzen (KNN)

ClearAll[f,interval,numberOfTrainingPoints,x,xTrain,yTrain,dataTrain,net,fKnn]

Generieren von Stützwerten anhand einer Beispielfunktion

f[x_]:=Exp[x Cos[3π x]]
interval={x,0,1};

Trainingsdatensätze erstellen

Im Allgemeinen Datensätze ist es vorteilhaft diese nach Trainings- und Testdaten aufzuteilen, um die Generalisierungsfähigkeit der KNN-Approximationen (Overfitting) überprüfen zu können

numberOfTrainingPoints=10;
xTrain=Range[interval[[2]],interval[[3]],((interval[[3]])/(numberOfTrainingPoints-1))]//N;
yTrain=Chop@(f[#]&/@xTrain);
dataTrain=Thread[xTrain->yTrain];

Definition der KNN-Struktur

(*net=NetChain[{5,ElementwiseLayer["Sigmoid"],5,ElementwiseLayer["Sigmoid"],1}]*)

Abkürzende Schreibweise:

net=NetChain[{5,Tanh,5,Tanh,1}]

Trainieren des KNN

fKnn=NetTrain[net,dataTrain]

Behandlung als “gewöhnliche” math. Funktion:

fKnn[(1/2)]

Plot[{f[x],fKnn[x]},interval,PlotStyle->{Blue,Green},PlotLegends->{"Funktion","Approx. KNN"},Epilog->Point[{xTrain,yTrain}\[Transpose]]]

Mit Wolfram Language Built-In Subroutinen

ClearAll[fPred]

Verschiedene Verfahren zur Interpolation sind in der Subroutine Predict zusammengefasst.

fPred=(Predict[dataTrain,Method->#]&)/@{"DecisionTree","GradientBoostedTrees","LinearRegression","NearestNeighbors","NeuralNetwork"}

Plot[{f[x],fPred[[1]][x],fPred[[2]][x],fPred[[3]][x],fPred[[4]][x],fPred[[5]][x]},interval,PlotStyle->{Blue,Green,Orange,Brown,Magenta,Red},PlotLegends->{"Funktion","Approx. Pred. D.-Tree","Approx. Pred. G. B.-Trees","Approx. Pred. Linear Regr.","Approx. Pred. Nearest Neighbor", "Approx. Pred. KNN"},PlotRange->All,Epilog->Point[{xTrain,yTrain}\[Transpose]]]

Daten Fitting und Generieren eines Modells mit verschieden Wolfram Language Built-In Subroutinen

Mit Vorgabe eines parametrisierten Modells

ClearAll[f,interval,numberOfPkts,xi,yi,pkts]

Generieren von Stützwerten anhand einer Beispielfunktion

f[x_]:=x Exp[2x]+6
interval={x,0,1};

numberOfPkts=10;
xi=Range[interval[[2]],interval[[3]],((interval[[3]])/(numberOfPkts-1))]//N;
yi=Chop@(f[#]&/@xi);
pkts={xi,yi}\[Transpose];

FindFit-Funktion

ClearAll[expr,a,b,c,sol]

Vorgabe einer parametrisierten Funktion

expr=a x Exp[b x]+c;

Fitten der Parameter auf vorgegebene Punkte

sol=FindFit[pkts,expr,{a,b,c},x,MaxIterations->1000,Method->NMinimize]

expr/. sol

Plot[{f[x],(expr/. sol)},interval,PlotStyle->{Blue,Green},PlotRange->All,Epilog->Point[pkts]]

NonlinearModelFit-Funktion

ClearAll[nlFittedModel]

nlFittedModel=NonlinearModelFit[pkts,expr,{a,b,c},x,MaxIterations->1000,Method->NMinimize]
Plot[nlFittedModel//Normal,interval,PlotStyle->Green,Epilog->Point[pkts]]

Zur Analyse der Genauigkeit und weiterer Daten-Fitting Modelle siehe Dokumentation Abschnitt Statistical Model Analysis

Ohne konkrete Vorgabe eines parametrisierten Modells

FindFormula[pkts,x]
Plot[%,interval,PlotStyle->Green,Epilog->Point[pkts]]

Vorgabe von zulässigen Basisfunktionen:

FindFormula[pkts,x,TargetFunctions->{Exp,Plus,Power}]
Plot[%,interval,PlotStyle->Green,Epilog->Point[pkts]]

Generieren verschiedener Lösungsvorschläge:

FindFormula[pkts,x,10,All,TargetFunctions->{Exp,Plus,Power}]

FindFormula[pkts,x,10,All]

Achtung: 
Vorhersagen nicht zuverlässig genau bei allen Funktionen (Stützwertverteilungen)

ClearAll[f,interval,numberOfPkts,xi,yi,pkts,expr,a,b,c,x]

f[x_]:=Exp[x Cos[3π x]];
interval={x,0,1};

numberOfPkts=10;
xi=Range[interval[[2]],interval[[3]],((interval[[3]])/(numberOfPkts-1))]//N;yi=Chop@(f[#]&/@xi);
pkts={xi,yi}\[Transpose];

Fitting ungenau, trotz Vorgabe des zugrundeliegenden math. Modells:

expr=a (E^(b x Cos[c x]));
NonlinearModelFit[pkts,expr,{a,b,c},x,MaxIterations->1000,Method->NMinimize]
FindFit[pkts,expr,{a,b,c},x,MaxIterations->1000,Method->NMinimize]
Plot[expr/.%,interval,Epilog->Point[pkts]]

Nur sehr grobe Annäherung an den vorgegebenen Datensatz:

FindFormula[pkts,x]
FindFormula[pkts,x,TargetFunctions->{Exp,Cos,Plus,Times,Power}]

Built-In Subroutine welche die verallgemeinerte math. Beschreibung von diskreten Folgengliedern zurückgibt

ClearAll[n]

FindSequenceFunction[{3,5,7,9,11,13,15,17,19,21},n]
FindSequenceFunction[{0,3,10,21,36,55,78,105,136,171},n]
FindSequenceFunction[{1,1,2,3,5,8,13},n]

FindGeneratingFunction[{0,3,10,21,36,55,78,105,136,171},x]