LAGRANGE Interpolation
Univ.-Prof. Dr.-Ing. habil. Josef BETTEN
RWTH Aachen University
Mathematical Models in Materials Science and Continuum Mechanics
Augustinerbach 4-20
D-52056 A a c h e n , Germany
<betten@mmw.rwth-aachen.de>
Durch n + 1 Stützwerte f (x[k]), k = 0,1,...,n, soll ein Polynom n-ten Grades gelegt werden. Mit
Hilfe der LAGRANGEschen Interpolationsmethode erhält man ein Polynom n-ten Grades gemäß:
restart:
P(x):=Sum(f(x[k])*L[n,k](x),k=0..n);
Für jedes k = 0,1, ..., n ist eine LAGRANGEsche Grundfunktion folgendermaßen definiert:
L[n,k](x):=Product((x-x[i])/(x[k]-x[i]),i=0..n), k<>i;
L[n,k](x[i]):=delta[ik]=piecewise(i=k,1, i<>k,0);
Beispiel: Gegeben seien sechs Messpunkte
[x[k],y[k]]=[[-1,0],[-1/2,3/4],[0,1],[1/4,10/9],[3/4,1/2],[1,0]];
durch die ein Interpolationspolynom fünften Grades gelegt werden soll.
Zunächst werden LAGRANGEsche Grundfunktionen fünften Grades (n = 5) formuliert.
L[5,k](x):=Product((x-x[i])/(x[k]-x[i]),i=0..5); k<>i;
L[5,k](x):=value(%%);
Für die einzelnen Stützstellen (k = 0,1,...,5) erhält man:
L[5,0](x):=%*(x[k]-x[0])/(x-x[0]);
L[5,0](x):= subs({x[k]=-1,x[1]=-1/2,x[2]=0,x[3]=1/4,x[4]=3/4,x[5]=1},%);
L[5,0](-1):=subs(x=-1,%);
L[5,0](x):=expand(%%);
L[5,1](x):=L[5,k](x)*(x[k]-x[1])/(x-x[1]);
L[5,1](x):= subs({x[k]=-1/2,x[0]=-1,x[2]=0,x[3]=1/4,x[4]=3/4,x[5]=1},%);
L[5,1](-1/2):=subs(x=-1/2,%);
L[5,1](x):=expand(%%);
L[5,2](x):=L[5,k](x)*(x[k]-x[2])/(x-x[2]);
L[5,2](x):= subs({x[k]=0,x[0]=-1,x[1]=-1/2,x[3]=1/4,x[4]=3/4,x[5]=1},%);
L[5,2](0):=subs(x=0,%);
L[5,2](x):=expand(%%);
L[5,3](x):=L[5,k](x)*(x[k]-x[3])/(x-x[3]);
L[5,3](x):= subs({x[k]=1/4,x[0]=-1,x[1]=-1/2,x[2]=0,x[4]=3/4,x[5]=1},%);
L[5,3](1/4):=subs(x=1/4,%);
L[5,3](x):=expand(%%);
L[5,4](x):=L[5,k](x)*(x[k]-x[4])/(x-x[4]);
L[5,4](x):= subs({x[k]=3/4,x[0]=-1,x[1]=-1/2,x[2]=0,x[3]=1/4,x[5]=1},%);
L[5,4](3/4):=subs(x=3/4,%);
L[5,4](x):=expand(%%);
L[5,5](x):=L[5,k](x)*(x[k]-x[5])/(x-x[5]);
L[5,5](x):= subs({x[k]=1,x[0]=-1,x[1]=-1/2,x[2]=0,x[3]=1/4,x[4]=3/4},%);
L[5,5](1):=subs(x=1,%);
L[5,5](x):=expand(%%);
Obige LAGRANGEsche Grundfunktionen können eleganter vermöge einer "DO-Schleife" direkt
erzeugt werden:
for i from 0 to 5 do L[5,i](x):=L[5,k](x)*(x[k]-x[i])/(x-x[i]) od;
Mit den gegebenen Stützwerten x[0],...,x[5] erhält man durch eine weitere "DO-Schleife"
die LAGRANGE Polynome:
for i from 0 to 5 do L[5,i](x):= subs({x[k]=x[i]}, {x[0]=-1,x[1]=-1/2,x[2]=0,x[3]=1/4,x[4]=3/4,x[5]=1},L[5,i](x)) od;
Weiterhin findet man:
for i from 0 to 5 do L[5,i](x):=expand(L[5,i](x)) od;
Im folgenden Bild sind die LAGRANGE Polynome L[5,k](x), k =0,1,...,5 dargestellt.
alias(th=thickness,co=color):
plot({seq(L[5,k](x),k=0..5)},x=-1..1,xtickmarks=4,th=3,co=black, title="LAGRANGE Funktionen L[5,k](x), k = 0. . . 5");
Mit obigen LAGRANGE Polynomen und den gegebenen Stützwerten f(x[k]) = y[k] ermittelt man
das Intepolationspolynom fünften Grades gemäß:
P[5](x):=Sum(y[kappa]*L[5,kappa](x),kappa=0..5);
with(linalg):
y[k]:=matrix(1,6,[0,3/4,1,10/9,1/2,0]);
L[5,k](x):=matrix(6,1,[seq(L[5,k](x),k=0..5)]);
P[5](x):=multiply(y[k],L[5,k](x));
P[5](x):=1+(4736/2835)*x^5-(1616/2835)*x^4-(6728/2835)*x^3- (1219/2835)*x^2+(664/945)*x;
Dieses Polynom mit den Messwerten
[x[k],y[k]]= [[-1,0],[-1/2,3/4],[0,1],[1/4,10/9],[3/4,1/2],[1,0]];
wird im nächsten Bild dargestellt.
alias(sc=scaling,th=thickness,co=color):
p[1]:=plot(P[5](x),x=-1..1,sc=constrained,th=3,co=black, title="Polynom fünften Grades"):
p[2]:=plot(rhs(%%%),x=-1..1,xtickmarks=4,ytickmarks=3, style=point,symbol=circle,symbolsize=30,th=4,co=black):
plots[display](seq(p[i],i=1..2));
Interpolationspolynome können auch mit Hilfe folgender Determinante erzeugt werden.
restart: with(linalg):
M:=Matrix([[P(x),seq(x^k,k=0..3),___,x^n], [y[0],seq(x[0]^k,k=0..3),___,x[0]^n], [y[1],seq(x[1]^k,k=0..3),___,x[1]^n], [y[2],seq(x[2]^k,k=0..3),___,x[2]^n], [°,°,°,°,°,°,°], [y[n],seq(x[n]^k,k=0..3),___,x[n]^n]]);
Das Polynom P(x) n-ten Grades erhält man aus der Determinante von M gemäß:
Polynom_P(x):=Solve(Det(M)=0,P(x));
Als Beispiel werde ein Interpolationspolynom fünften Grades gesucht durch die oben
angenommenen sechs Knotenpunkte.
M:=matrix(7,7,[[P[5](x),seq(x^k,k=0..5)], [y[0],seq(x[0]^k,k=0..5)], [y[1],seq(x[1]^k,k=0..5)], [y[2],seq(x[2]^k,k=0..5)], [y[3],seq(x[3]^k,k=0..5)], [y[4],seq(x[4]^k,k=0..5)], [y[5],seq(x[5]^k,k=0..5)]]);
data:=[-1,0],[-1/2,3/4],[0,1],[1/4,10/9],[3/4,1/2],[1,0];
m:=subs({x[0]=-1,x[1]=-1/2,x[2]=0,x[3]=1/4,x[4]=3/4,x[5]=1}, {y[0]=0,y[1]=3/4,y[2]=1,y[3]=10/9,y[4]=1/2,y[5]=0},%%);
P[5](x):=solve(det(m)=0,P[5](x));
Dieses Ergebnis stimmt mit der LAGRANGE Interpolation überein.
Mit Maple erhält man durch den Befehl evalf eine Gleitpunktdarstellung (floating point arithmetic) gemäß.
Q[5](x):=evalf(P[5](x));
Delta(x):=P[5](x)-Q[5](x);
Maple rechnet mit einer Mantissenlänge von 10. Man erhält im obigen Beispiel das exakte Ergebnis.
Je nach Problem kann die Genauigkeit erhöht werden durch eine zusätzliche Angabe: Digits > 10.
Ebenso erhält man auch Näherungen durch Digits < 10, wie folgendermaßen gezeigt wird:
q[5](x):=evalf(P[5](x),3);
Als Abweichung gilt:
delta(x):=P[5](x)-q[5](x);
plot(delta(x),x=-1..1,th=3,co=black,xtickmarks=4,ytickmarks=3);
Fehlernorm
L[2]:=sqrt((1/2)*Int((delta(xi))^2,xi=-1..1))= evalf(sqrt((1/2)*int((delta(x))^2,x=-1..1)),4);
Bei einer Vielzahl von Knotenpunkten (Meßpunkten) führt die LAGRANGE Interpolation auf ein
Polynom hohen Grades, das je nach Lage der Stützstellen zu mehr oder weniger starken "Schwingungen" neigt. Die LAGRANGE Interpolation ist in solchen Fällen nicht sehr nützlich. Vorteilhafter sind dann interpolierende Splines zur Erzeugung glatter Kurven, wie im Folgenden gezeigt wird.
restart: with(CurveFitting):
S[3](x):=Spline([data],x,degree=3);
P[5](x):=1+(664/945)*x-(1219/2835)*x^2-(6728/2835)*x^3- (1616/2835)*x^4+(4736/2835)*x^5;
p[1]:=plot(S[3](x),x=-1..1,th=3,co=black, title="Polynom fünften Grades und kubischer Spline"):
p[2]:=plot(P[5](x),x=-1..1,ytickmarks=3,th=2,co=black):
p[3]:=plot([data],style=point,symbol=circle,symbolsize=30):
plots[display](seq(p[k],k=1..3));
Man erkennt die "schwingende" Form des LAGRANGE Polynoms fünften Grades, während die Spline Interpolation stückweise kubische Polynome erzeugt.
Die Bezeichnung "Spline" bedeutet "dünnes Brett". Ein leicht biegsames Lineal wird an die
Knotenpunkte gelegt. Somit kann eine geglättete Kurve gezeichnet werden.
Im obigen Bild sind ein Interpolationspolynom P[5](x) fünften Grades und eine kubische
Splinefunktion dargestellt. Ergänzend sind im folgenden Bild ein linearer Spline und eine Splinefunktion fünften Grades ausgedruckt.
for i in [1,5] do S[i](x):=Spline([data],x,degree=i) od;
p[1]:=plot(S[1](x),x=-1..1,th=1,co=black, title="P[5](x), S[1](x), S[3](x), S[5](x)"):
p[2]:=plot(S[3](x),x=-1..1,ytickmarks=3,th=2,co=black):
p[3]:=plot(S[5](x),x=-1..1,th=1,co=black):
p[4]:=plot(P[5](x),x=-1..1,th=3,co=black):
p[5]:=plot([data],style=point,symbol=circle, symbolsize=30,th=3,co=black):
plots[display](seq(p[k],k=1..5));
Darin verbindet der lineare Spline die Knotenpunkte linear. Die Splinefunktionen dritten und fünften Grades unterscheiden sich geringfügig. Sie bilden eine gute Intepolation der Messpunkte. Wegen des geringeren Aufwandes ist die kubische Splineinterpolation vorzuziehen. Darüber
hinaus schmiegen sich die Spline-Kurven mit wachsendem " degree" immer stärker an das
"schwingende" Interpolations Polynom P[5](x). Im Grenzfall degree ---> infinity sind das
LAGRANGEsche Interpolations Polynom und die Spline-Kurve deckungsgleich:
with(CurveFitting):
for i in [1,3,5,10,11] do S[i](x):=Spline([data],x,degree=i) od:
p[1]:=plot({S[1](x),S[5](x),S[10](x),S[11](x)},x=-1..1, sc=constrained,xtickmarks=4,ytickmarks=3,th=1,co=black):
p[2]:=plot(S[3](x),x=-1..1,th=2,co=black):
p[3]:=plot(P[5](x),x=-1..1,th=3,co=black, title="P[5], S[1], S[3], S[5], S[10], S[11]"):
p[4]:=plot([data],style=point,symbol=circle,symbolsize=30, th=4,co=black):
plots[display](seq(p[k],k=1..4));
Abstand zwischen P[5](x) und S[n](x) ausgedrückt durch die Fehlernorm L-zwei
L[2]:=sqrt((1/2)*Int((P[5](xi)-S[n](xi))^2,xi=-1..1));
for i in [1,3,5,10,11] do L[2][i]:=evalf(sqrt((1/2)*int((P[5](x)-S[i](x))^2,x=-1..1))) od;
Grenzwert für dgree = infinity
L[2][infinity]:= Limit(sqrt((1/2)*Int((P[5](xi)-S[n](xi))^2,x=-1..1)),n=infinity)=0;
Bereits für degree >= 11 fällt die Splinefunktion S[n](x) im hier betrachteten Beispiel mit dem
LAGRANGEschen Interpolationspolynom P[5](x) nahezu zusammen.
Formfunktionen (shape functions) für finite Elemente
Die LAGRANGEsche Interpolationsmethode ist neben anderen Verfahren sehr hilfreich zur
Erzeugung von Formfunktionen sowohl für ebene als auch für räumliche finite Elemente.
Formfunktionen N(x,y,z) sind Interpolationsfunktionen, die einen Zusammenhang zwischen den
Knotenvariablen u[k], ... , w[k] und den Verschiebungen u(x,y,z), ..., w(x,y,z) im Innern eines
finiten Elementes herstellen gemäß:
[u(x,y,z),__,w(x,y,z)]= [Sum(u[k]*N[k](x,y,z),k=1..n),__,Sum(w[k]*N[k](x,y,z),k=1..n)];
Beispielsweise werden in [BETTEN, J.: Finite Elemente für Ingenieure 1, zweite Aufl., 2003,
Springer, Berlin / Heidelberg] eine Vielzahl von Formfunktionen der LAGRANGE - Klasse
und anderer Typen für Dreiecks-, Rechteck-, Tetraeder-, Hexaederelemente etc. aufgestellt.
