Interpolation - Beuth Hochschule für Technik...
Transcript of Interpolation - Beuth Hochschule für Technik...
Interpolation Ausgehend von N Datenpunkten (yk, xk) = {yk, xk | k = 1, 2, …, N} (Bild 1, a) wird bei der Interpo-lation eine mathematische Funktion y = f(x) gesucht, die einerseits die gegebenen Datenpunkte miteinander verbindet und andererseits eine einfache mathematische Beziehung liefert mit der man auch Zwischenwerte berechnen kann. Hierfür stellt die numerische Mathematik eine Vielzahl von Ansätzen und Verfahren zur Verfügung, von denen im Folgenden die lineare Interpolation sowie die Interpolation auf der Basis der kubischen Splines behandelt werden. 1 Lineare Interpolation Die lineare Interpolation ist ein sehr einfaches Verfahren und wird deshalb in der Praxis sehr häufig benutzt. Hierbei werden die einzelnen Datenpunkte durch Geraden miteinander verbunden (Bild 1, b). Ausgangspunkt ist die allgemeine Geradengleichung: y = y0 + m ⋅ x (y0: Achsabschnitt, m: Steigung). Sind nun N Datenpunkte (yk, xk) = {yk, xk | k = 1, 2, …, N} gegeben, so lassen sich N - 1 Geraden-gleichungen aufstellen:
1-Nk1für ; )xx(xxyyyy k
k1k
k1kk ≤≤−⋅
−−
+=+
+ und xk < x < xk+1
Lineare Interpolation mit MATLAB
Die graphische Darstellung hierzu zeigt Bild 1. b.
% Lineare Interpolation x = [-0.5 1 1.5 2.3 2.9 4.1]; % x-Datenwerte y = [1.2 0.25 -0.8 -2.3 0.8 -0.26]; % y-Datenwerte xi = -0.5:.1:4.1; % x-Interpolationswerte yi = interp1(x,y,xi,'linear'); % y-Interploationswerte plot(x,y,'o',xi,yi),axis([-1,4.5,-3.5,4])
-1 0 1 2 3 4
-2
0
2
4
-1 0 1 2 3 4
-2
0
2
4
-1 0 1 2 3 4
-2
0
2
4
Bild 1. a) Datenpunkte b) lineare Interpolation c) kubische Spline Interpolation
a) b) c)
2
2 Interpolation mit kubischen Splines Die lineare Interpolation ergibt zwar insgesamt einen stetigen, aber keinen “glatten“ Verlauf der Inter-polationsfunktion. Glatt bedeutet dabei, dass die zweite Ableitung der Interpolationsfunktion eine kontinuierliche Funktion ist. Weiterhin soll die Funktion außerhalb des Intervalls [x1, xN] ein bestimm-tes Verhalten (Randbedingungen) aufweisen. Derartige Eigenschaften hat eine Interpolationsfunktion auf der Basis der kubischen Splines (Bild 1.c). Setzt man wiederum N Datenpunkte (yk, xk) = {yk, xk | k = 1, 2, …, N} voraus (im Folgenden auch Stützstellen genannt), so wird hierbei für jedes Teilinter-vall [x1, x2], [x2, x3], …, [xN-1, xN] ein Polynom 3. Grades bestimmt, womit sich insgesamt N-1 derar-tige Polynome ergeben:
3kk
2kkkkkk )xx(d)xx(c)xx(ba)x(P −⋅+−⋅+−⋅+= ; (1)
für xk < x < xk+1 und 1 < k < N-1.
Hierbei wird durch die Bedingung (a) sichergestellt, dass die Splinefunktion alle Stützstellenpunkte miteinander verbindet und durch die Bedingungen (b, c) ist gewährleistet, dass Interpolationsfunktion insgesamt einen glatten Verlauf zeigt. Sind N Datenpunkte gegeben so sind N-1 Polynome mit jeweils 4 Parametern, insgesamt also 4⋅(N-1) unbekannte Parameter zu bestimmen. Zur Bestimmung dieser Unbekannten werden die Bedingungen (a – d) herangezogen. Dabei ergibt sich aus der Bedingung (a) für die Parameter ak die Beziehung: ak = yk ; für 1 < k < N-1 (2) Für die übrigen Parameter haben die Parameter ck der quadratischen Terme eine besondere Bedeutung, da sich die Parameter bk und dk hiervon ableiten lassen. Eine ausführliche mathematische Herleitung dieser Zusammenhänge ist dem Anhang zu entnehmen. Im Folgenden werden die dabei erzielten Ergebnisse gezeigt und in Verbindung mit einem Beispiel angewendet. Definiert man als:
Parametervektor:
⎟⎟⎟⎟⎟⎟⎟⎟⎟
⎠
⎞
⎜⎜⎜⎜⎜⎜⎜⎜⎜
⎝
⎛
⋅⋅⋅
=
−
N
1N
2
1
cc
cc
Cr
(3)
Jedes der Polynome muß an den Stützstellen die folgenden Bedingungen erfüllen: a) Pk(xk) = yk und Pk(xk+1) = yk+1 b) 1. Ableitung: )x(P)x(P 1k
'1k1k
'k +++ =
c) 2. Ableitung: )x(P)x(P 1k
''1k1k
''k +++ =
d) Randbedingung: 0)x(P)x(P N
''1N1
''1 == − (natürliche Splines)
3
Hierbei ist cN = 0 (Pseudoparameter). Mit den Differenzwerten: k1kk xx −=ε + und k1kk yy −=δ + definiert man weiterhin den
Differenzenvektor:
⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟
⎠
⎞
⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜
⎝
⎛
εδ
−εδ
εδ
−εδ
εδ
−εδ
⋅=
−
−
−
−
0
:
0
3Y
2N
2N
1N
1N
2
2
3
3
1
1
2
2
r (4)
und die −× NN Differenzenmatrix: Damit läßt sich die folgende Matrizengleichung aufstellen:
YCrr
=⋅M (6) Bildet man die inverse Differenzenmatrix M-1, so ergibt sich damit der Parametervektor gemäß:
YC 1rr
⋅= −M (7) Mit den N Elementen ck des Parametervektors C
rlassen sich dann die Polynomparameter bk gemäß:
k1kk31
k
kk )cc2(b ε⋅+⋅⋅−
εδ
= + ; für 1 < k < N-1 (8)
und dk gemäß:
k
k1kk 3
ccdε⋅−
= + ; für 1 < k < N-1; mit cN = 0 (9)
berechnen. 2. Beispiel: Gegeben sind die nebenstehenden 5 Datenpunkte (xk, yk), für die eine Interpolationsfunktion auf der Basis der natürlichen kubischen Splines aufzustel-len ist.
Datenpunkte
k xk yk εk = xk+1 - xk δk = yk+1 - yk 1 2 3 4 5
1 2,5 3,3 4
4,5
1 -0,5
2 -1 1,1
1,5 0,8 0,7 0,5
-1,5 2,5 -3 2,1
⎟⎟⎟⎟⎟⎟⎟⎟⎟
⎠
⎞
⎜⎜⎜⎜⎜⎜⎜⎜⎜
⎝
⎛
εε+ε⋅ε
εε+ε⋅εεε+ε⋅ε
=
−−−−
100...0000)(2...0000
000...)(20000...0)(2000...0001
M
1N1N2N2N
3322
2211
(5)
4
Hierzu sind 4 kubische Polynome:
3kk
2kkkkkk )xx(d)xx(c)xx(ba)x(P −⋅+−⋅+−⋅+= ; mit 1 < k < 4 und x ∈ [xk, xk+1]
erforderlich. Jedes weist 4 Parameter auf, womit insgesamt 16 unbekannte Parameter zu bestimmen sind. Nach (2) ergeben sich die vier Parameter ak zu: ak = yk ; für 1 < k < 4 Zur Bestimmung der Parameter ck; für 1 < k < 5 werden nach (3) der Parametervektor:
⎟⎟⎟⎟⎟⎟
⎠
⎞
⎜⎜⎜⎜⎜⎜
⎝
⎛
=
5
4
3
2
1
ccccc
Cr
; mit c5 = 0
nach (4) der Differenzenvektor:
⎟⎟⎟⎟⎟⎟
⎠
⎞
⎜⎜⎜⎜⎜⎜
⎝
⎛
−=
⎟⎟⎟⎟⎟⎟⎟⎟⎟
⎠
⎞
⎜⎜⎜⎜⎜⎜⎜⎜⎜
⎝
⎛
εδ
−εδ
εδ
−εδ
εδ
−εδ
⋅=
04571,252321.22375,120
0
0
3Y
3
3
4
4
2
2
3
3
1
1
2
2
r
und nach (5) die Differenzenmatrix:
⎟⎟⎟⎟⎟⎟
⎠
⎞
⎜⎜⎜⎜⎜⎜
⎝
⎛
=
⎟⎟⎟⎟⎟⎟
⎠
⎞
⎜⎜⎜⎜⎜⎜
⎝
⎛
εε+ε⋅εεε+ε⋅ε
εε+ε⋅ε=
100005,04,27,000
07,038,00008,06,45,100001
10000)(200
0)(2000)(200001
M
4343
3232
2121
aufgestellt. Aus M ergibt sich die inverse Differenzenmatrix:
⎟⎟⎟⎟⎟⎟
⎠
⎞
⎜⎜⎜⎜⎜⎜
⎝
⎛
−−−−
−−=−
100002243,04487,01098.00191,00286,0
0549,01098,03764,00655,00982,00095,00191,00655,02288,03432,000001
1M
womit sich dann nach (7) der Parametervektor:
5
⎟⎟⎟⎟⎟⎟
⎠
⎞
⎜⎜⎜⎜⎜⎜
⎝
⎛
−=
⎟⎟⎟⎟⎟⎟
⎠
⎞
⎜⎜⎜⎜⎜⎜
⎝
⎛
00993,149732,11
7725,40
ccccc
5
4
3
2
1
ergibt. Mit diesen Werten berechnen sich dann mit Gl. 8 und 9 die Parameter bk und dk, womit sich insgesamt die folgenden Polynomparameter:
k a b c d 1 2 3 4
1 -0,5
2 -1
-3,3863 3,7725 -1,9881 -0,4998
0 4,7725
-11,9732 14,0993
1,0606 -6,9774 12,4155 -9,3996
und damit die folgenden Polynome (Bild 4) ergeben:
Das Polynom P1(x) gilt auch für Werte x < 1 und P4(x) gilt auch für Werte x > 4,5 also für Werte außerhalb der Stützstellen (Extrapolation). Zu bemerken ist: Vergleicht man verschiedene Spline-Interpolationsverfahren, so können diese aufgrund von unterschiedlichen Randbedingungen in den Randbereichen und insbesondere darüber
für 1 < x < 2,5 : 31 )1x(061,1)1x(386,31)x(P −⋅+−⋅−=
für 2,5 < x < 3,3: 32
2 )5,2x(977,6)5,2x(773,4)5,2x(773,35,0)x(P −⋅−−⋅+−⋅+−= ; für 3,3 < x < 4: 32
3 )3,3x(416,12)3,3x(973,11)3,3x(989,12)x(P −⋅+−⋅−−⋅−= ; für 4 < x < 4,5: 32
4 )4x(4,9)4x(099,14)4x(5,01)x(P −⋅−−⋅+−⋅−−=
Bild 4. Datenpunkte und Spline-Funktion
0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5-1.5
-1
-0.5
0
0.5
1
1.5
2
2.5
3
3.5
x →
↑y
6
hinaus, d.h. für Werte x < x1 und x > xN , abweichende Ergebnisse liefern. So wird auch bei MATLAB eine gegenüber (d) andere Randbedingung benutzt. Näheres hierzu ist ebenfalls dem Anhang zu ent-nehmen. 3 Kubische Spline-Interpolation mit MATLAB Für eine Interpolation von Messwerten mit kubischen Splines hat MATLAB die Funktion spline implementiert. Der Einsatz dieser Funktion wird mit folgendem Beispiel gezeigt. 3. Beispiel Gegeben sind die nebenstehenden Datenpunkte (identisch mit denen des 2. Beispiels), für die eine Interpolationsfunktion auf der Basis der kubischen Splines graphisch darzustellen ist. MATLAB-Programm zur Interpolation der Datenpunkte durch kubische Splines:
Bild 5 zeigt die Interpolationsfunktion mit den fünf Datenpunkten. Gegenüber der Randbedingung (d) in Kapitel 2 benutzt MATLAB eine andere Randbedingung, wobei eine Kontinuität der 3. Ableitung an der zweiten und vorletzten Stützstelle gefordert wird. Im vorlie-genden Beispiel bedeutet das für die Stützstelle x2 : )x(P)x(P 2
'''22
'''1 = und für die Stützstelle x4:
0.5 1 1.5 2 2.5 3 3.5 4 4.5 5-5
-4
-3
-2
-1
0
1
2
3
4
5
x →
↑y
Bild 5. MATLAB-Plot
k xk yk 1 2 3 4 5
1 2,5 3,3 4
4,5
1 -0,5
2 -1 1,1
Datenpunkte
% Interpolation mit kubischen Splinesx = [1 2.5 3.3 4 4.5]; % x-Datenwerte y = [1 -0.5 2 -1 1.1]; % y-Datenwerte xi = 0:.01:5; % x-Interpolationswerte yi = spline(x,y,xi); % y-Interpolationswerte plot(x,y,'ro',xi,yi,'b'), axis([0.5,5,-5,5])
7
)x(P)x(P 4'''
44'''
3 = . Hierdurch treten in der Umgebung der Randpunkte zwischen den beiden Splines erhebliche Abweichungen auf, wie Bild 6 zeigt. Für Interessierte werden die mathematischen Grundlagen hierzu im Anhang gezeigt. Anhang Gegeben sind N Datenpunkte (xk, yk); für 1 < k < N, auch Stützstellen genannt, aus denen sich für jeweils zwei unmittelbar aufeinander folgende Werte die Differenzen der x-Werte εk und der y-Werte δk berechnen lassen: Zur Interpolation der N Datenpunkte sind N-1 kubische Polynome:
0.5 1 1.5 2 2.5 3 3.5 4 4.5 5-5
-4
-3
-2
-1
0
1
2
3
4
5
x →
↑y
MATLAB-Splines
Bild 6. Vergleich verschiedener Splines
natürlicheSplines
k xk yk εk = xk+1 - xk δk = yk+1 - yk 1 2 :
N-1 N
x1 x2 :
xN-1 xN
y1 y2 :
yN-1 yN
ε1 = x2 – x1 ε2 = x3 – x2
: εN-1 = xN - xN-1
δ1 = y2 – y1 δ2 = y3 – y2
: δN-1 = yN - yN-1
x1 x2 x3 xN-2 xN-1 xN
P1
P2
PN-2
PN-1
8
3kk
2kkkkkk )xx(d)xx(c)xx(ba)x(P −⋅+−⋅+−⋅+= ; (A1)
mit 1 < k < N-1 und x ∈ [xk, xk+1] erforderlich. Jedes Polynom hat 4 Parameter, womit insgesamt 4⋅(N - 1) Parameter zu bestimmen sind. Hierzu werden die folgenden drei allgemeinen Bedingungen (a, b, c) herangezogen, mit denen sich dann 4⋅N-6 Gleichungen aufstellen lassen. Allgemeine Bedingungen: a) Pk(xk) = yk und Pk(xk+1) = yk+1 b) 1. Ableitung: )x(P)x(P 1k
'1k1k
'k +++ =
c) 2. Ableitung: )x(P)x(P 1k
''1k1k
''k +++ =
Zur Aufstellung der noch fehlenden 2 Gleichungen gibt es verschiedene Möglichkeiten, beispielsweise in eine der folgenden Randbedingungen: d1) natürliche Splines: 0)x(P)x(P N
''1N1
''1 == −
d2) Kontinuität der 3. Ableitung an der zweiten und vorletzten Stützstelle („not-a-knot“ Spline, implementiert bei der spline-Funktion von MATLAB): )x(P)x(P 2
'''22
'''1 = und )x(P)x(P 2N
'''2N2N
'''2N −−−− =
d3) Splines mit eingespanntem Rand: α=)x(P 1
'1 und β=− )x(P N
'1N
d4) Splines mit periodischer Randbedingung: )x(P)x(P N
'1N1
'1 −= und )x(P)x(P N
''1N1
''1 −=
Aufgrund der gewählten Randbedingung können die einzelnen Spline-Funktionen, insbesondere in den Randbereichen, erheblich voneinander abweichen. Die ersten drei Ableitungen der kubischen Polynome betragen:
2kkkkk
'k )xx(d3)xx(c2b)x(P −⋅⋅+−⋅⋅+= (A2)
)xx(d6c2)x(P kkk
''k −⋅⋅+⋅= (A3)
k'''
k d6)x(P ⋅= (A4) (i) Allgemeine Bedingungen Berücksichtigt man zunächst nur die drei allgemeinen Bedingungen, so ergeben sich folgende Bestim-mungsgleichungen: Aus der Bedingung (a) und Gl. A1 erhält man:
kkkk ya)x(P == Damit gelten für die Polynomparameter: ak = yk ; für 1 < k < N-1 (A5)
9
Weiterhin ergibt sich aus (a):
1k3
k1kk2
k1kkk1kkk1kk y)xx(d)xx(c)xx(ba)x(P +++++ =−⋅+−⋅+−⋅+= (A6) Berücksichtigt man dabei (A5), so folgt:
1k3
k1kk2
k1kkk1kkk y)xx(d)xx(c)xx(by ++++ =−⋅+−⋅+−⋅+ ; für 1 < k < N woraus sich nach einer elementaren Umformung:
k1k
k1k2k1kkk1kkk xx
yy)xx(d)xx(cb−−
=−⋅+−⋅++
+++
ergibt. Definiert man zur Abkürzung der Schreibweise die beiden Differenzwerte:
k1kk xx −=ε + und
k1kk yy −=δ + so erhält man die Beziehung:
k
k2kkkkk dcb
εδ
=ε⋅+ε⋅+ ; für 1 < k < N-1 (A7)
Mit Gl. A2 gilt:
2k1kkk1kkk1k
'k )xx(d3)xx(c2b)x(P −⋅⋅+−⋅⋅+= +++
und
1k1k'
1k b)x(P +++ = Berücksichtigt man die Bedingung (b), so erhält man:
2kkkkk1k d3c2bb ε⋅⋅+ε⋅⋅+=+ ; für 1 < k < N-2 (A8)
Weiterhin gilt mit Gl. A3:
)xx(d6c2)x(P k1kkk1k''
k −⋅⋅+⋅= ++ und
1k1k''
1k c2)x(P +++ ⋅= Damit folgt aus der Bedingung (c) die Beziehung:
kkk1k d3cc ε⋅⋅+=+ ; für 1 < k < N-2 (A9) Resümee: Mit den Beziehungen (A5, A7, A8 und A9) lassen sich insgesamt 4⋅N – 6 Gleichungen auf-stellen. Da jedoch 4⋅N – 4 Parameter zu bestimmen sind, fehlen noch 2 Gleichungen. Diese ergeben sich jeweils aus den folgenden Randbedingungen.
10
(ii) Splines mit der Randbedingung d1 (natürliche Splines) Berücksichtigt man die Randbedingung für natürliche Splines sowie die Gl. A3, so ergeben sich die beiden Beziehungen:
0c)x(P 11''
1 == (A10) und
0d6c2)xx(d6c2)x(P 1N1N1N1NN1N1NN''
1N =ε⋅⋅+⋅=−⋅⋅+⋅= −−−−−−− Aus der letzteren folgt:
0d3c 1N1N1N =ε⋅⋅+ −−− Für den weiteren Rechengang ist es zweckmäßig, mit dieser Beziehung neben den N-1 regulären Polynomparametern c1, c2, …, cN-1 einen weiteren Parameter, den Pseudoparameter cN einzuführen, gemäß:
0cd3c N1N1N1N ==ε⋅⋅+ −−− (A11) Damit erweitert sich dann der Geltungsbereich von Gl. A9 um den Wert k = N-1 und lautet nunmehr: kkk1k d3cc ε⋅⋅+=+ ; für 1 < k < N-1 Stellt man diese Gleichung nach dk um, so ergibt sich:
k
k1kk 3
ccdε⋅−
= + ; für 1 < k < N-1; mit cN = 0 (A12)
Damit lassen sich die N-1 Polynomparameter dk mit den N Parametern c1, c2, …, cN-1, cN berechnen, wobei c1 = cN = 0. Im weiteren wird nun gezeigt, dass sich auch die N-1 Polynomparameter bk mit den N Werten ck berechnen lassen. Berücksichtigt man in (A6) für dk die Beziehung (A12), so ergibt sich:
k
kk1kk3
1k )cc2(b
εδ
=ε⋅+⋅⋅+ + (A13)
Hieraus folgt:
k1kk31
k
kk )cc2(b ε⋅+⋅⋅−
εδ
= + ; für 1 < k < N-1 (A14)
Damit lassen sich auch die N-1 Parameter bk mit den N Parametern c1, c2, …, cN-1, cN berechnen. Insofern haben bei den natürlichen Splines die Parameter ck eine zentrale Bedeutung. Nun zur Bestimmung der Parameter ck; für 1 < k < N. Hierzu müssen N Bestimmungsgleichungen aufgestellt werden. Ausgang hierzu ist Gl. A7. Berücksichtigt man in (A7) für dk (A12), so erhält man:
kk1kk1k )cc(bb ε⋅++= ++ ; für 1 < k < N-2
11
Nimmt man hierbei für den Index k die Substitution: k → k – 1 vor, so gilt:
1k1kk1kk )cc(bb −−− ε⋅++= ; für 2 < k < N-1 Berücksichtigt man hierbei für die beiden Parameter bk und bk-1 (A14), so ergibt sich:
1k1kk1kk1k31
1k
1kk1kk3
1
k
k )cc()cc2()cc2( −−−−−
−+ ε⋅++ε⋅+⋅⋅−
εδ
=ε⋅+⋅⋅−εδ
; für 2 < k < N-1
woraus sich nach einigen elementaren Umformungen N-2 Bestimmungsgleichungen für die N Para-meter ck ergeben, der Form:
3 c)(c2c1k
1k
k
kk1k1kkk1k1k ⎟⎟
⎠
⎞⎜⎜⎝
⎛εδ
−εδ
⋅=ε⋅+ε+ε⋅⋅+ε⋅−
−+−−− ; für 2 < k < N-1
Berücksichtigt man noch die beiden Beziehungen c1 = cN = 0 nach (A9) und (A10), so ergeben sich zwei weitere Gleichungen, womit sich schließlich das folgende System mit insgesamt N Bestim-mungsgleichungen für die Parameter ck ergibt:
c1 = 0
3 c)(c2c1k
1k
k
kk1k1kkk1k1k ⎟⎟
⎠
⎞⎜⎜⎝
⎛εδ
−εδ
⋅=ε⋅+ε+ε⋅⋅+ε⋅−
−+−−− ; 2 < k < N-1
cN = 0 Dieses Gleichungssystem läßt sich durch Einführung der folgenden Definitionen: (i) Parametervektor:
⎟⎟⎟⎟⎟⎟⎟⎟⎟
⎠
⎞
⎜⎜⎜⎜⎜⎜⎜⎜⎜
⎝
⎛
⋅⋅⋅
=
−
N
1N
2
1
cc
cc
Cr
(ii) Differenzenvektor:
⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟
⎠
⎞
⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜
⎝
⎛
εδ
−εδ
εδ
−εδ
εδ
−εδ
⋅=
−
−
−
−
0
:
0
3Y
2N
2N
1N
1N
2
2
3
3
1
1
2
2
r
12
(iii) Differenzenmatrix: als Matrizengleichung formulieren:
YCrr
=⋅M Bildet man die inverse Differenzenmatrix M-1, so ergibt sich damit der Parametervektor gemäß:
YC 1rr
⋅= −M Mit den N Elementen des Parametervektors lassen sich dann mit (A14) die Polynomparameter bk und mit (A12) die Parameter dk berechnen. Berücksichtigt man noch Gl. A4, so sind damit sämtliche Parameter der N-1 kubischen Polynome bekannt. (ii) Splines mit der Randbedingung d2 Bei dieser Randbedingung (bekannt auch unter der Bezeichnung „not-a-knot“ Spline) wird eine Kontinuität der 3. Ableitung an der zweiten und vorletzten Stützstelle gefordert. Zu bemerken ist, dass diese Randbedingung bei der spline-Funktion von MATLAB implementiert ist. Damit ergibt sich für die Stützstellen: x2: )x(P)x(P 2
'''22
'''1 =
xN-1: )x(P)x(P 1N
'''1N1N
'''2N −−−− = .
Aus (A3) folgt für die 3. Ableitung: k
'''k d6P ⋅= ; für 1 < k < N-1. Damit ergeben sich mit der
Randbedingung (d2) die folgenden Beziehungen:
)x(P)x(P 2'''
22'''
1 = ⇒ d1 = d2 (A15) und
)x(P)x(P 1N'''
1N1N'''
2N −−−− = ⇒ dN-2 = dN-1 (A16) Mit d1 = d2 folgt aus (A12):
2
23
1
12
3cc
3cc
ε⋅−
=ε⋅−
woraus sich die Beziehung:
⎟⎟⎟⎟⎟⎟⎟⎟⎟
⎠
⎞
⎜⎜⎜⎜⎜⎜⎜⎜⎜
⎝
⎛
εε+ε⋅ε
εε+ε⋅εεε+ε⋅ε
=
−−−−
100...0000)(2...0000
000...)(20000...0)(2000...0001
M
1N1N2N2N
3322
2211
13
0cc)(c 3122112 =⋅ε+⋅ε+ε−⋅ε (A17) ergibt. Mit dN-2 = dN-1 folgt aus (A1):
N3
1NN2N2
1NN1N1NN1N1N y)xx(d)xx(c)xx(ba =−⋅+−⋅+−⋅+ −−−−−−− Hieraus ergibt sich unter Berücksichtigung von (A5):
1N
1N21N2N1N1N1N dcb
−
−−−−−− ε
δ=ε⋅+ε⋅+ (A18)
Setzt man in (A7) für k = N-2, so erhält man:
22N2N2N2N2N1N d3c2bb −−−−−− ε⋅⋅+ε⋅⋅+=
Dies eingesetzt in (A18) ergibt:
( )1N
1N2N
21N
22N1N1N2N2N2N d3cc2b
−
−−−−−−−−− ε
δ=⋅ε+ε⋅+ε⋅+ε⋅⋅+
Setzt man in (A12) für k = N-2, so erhält man:
2N
2N1N2N 3
ccd−
−−− ε⋅
−=
Dies eingesetzt in die vorangehende Gleichung ergibt:
( )1N
1N
2N
2N1N21N
22N1N1N2N2N2N 3
cc3cc2b−
−
−
−−−−−−−−− ε
δ=
ε⋅−
⋅ε+ε⋅+⋅ε+⋅ε⋅+
Setzt man in (A14) für k = N-2 so erhält man:
2N1N2N31
2N
2N2N )cc2(b −−−
−
−− ε⋅+⋅⋅−
εδ
=
Setzt man dies in die vorangehende Gleichung ein, so erhält man nach Umformungen:
⎟⎟⎠
⎞⎜⎜⎝
⎛εδ
−εδ
⋅=⎟⎟⎠
⎞⎜⎜⎝
⎛εε
+ε⋅+ε⋅+⋅⎟⎟⎠
⎞⎜⎜⎝
⎛εε
−ε−
−
−
−
−
−−−−
−
−−
2N
2N
1N
1N
2N
21N
2N1N2N2N
21N
2N 323c
Damit ergibt sich zur Berechnung der N-1 Parameter ck das Gleichungssystem:
0cc)(c 3122112 =⋅ε+⋅ε+ε−⋅ε
3 cc)(2c1k
1k
k
k1kkk1kk1k1k ⎟⎟
⎠
⎞⎜⎜⎝
⎛εδ
−εδ
⋅=⋅ε+⋅ε+ε⋅+⋅ε−
−+−−− ; für 2 < k < N-2
⎟⎟⎠
⎞⎜⎜⎝
⎛εδ
−εδ
⋅=⎟⎟⎠
⎞⎜⎜⎝
⎛εε
+ε⋅+ε⋅+⋅⎟⎟⎠
⎞⎜⎜⎝
⎛εε
−ε−
−
−
−
−
−−−−
−
−−
2N
2N
1N
1N
2N
21N
2N1N2N2N
21N
2N 323c
14
Führt man wiederum die Definitionen: (i) Parametervektor:
⎟⎟⎟⎟⎟⎟⎟⎟⎟
⎠
⎞
⎜⎜⎜⎜⎜⎜⎜⎜⎜
⎝
⎛
⋅⋅⋅
=
−
−
1N
2N
2
1
cc
cc
Cr
(ii) Differenzenvektor:
⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟
⎠
⎞
⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜
⎝
⎛
εδ
−εδ
εδ
−εδ
εδ
−εδ
εδ
−εδ
⋅=
−
−
−
−
−
−
−
−
2N
2N
1N
1N
3N
3N
2N
2N
2
2
3
3
1
1
2
2
:
0
3Yr
(iii) Differenzenmatrix:
ein, so läßt die Matrizengleichung:
YCrr
=⋅M formulieren aus der sich der Parametervektor gemäß:
YC 1rr
⋅= −M ergibt. Die Parameter ak ergeben sich wiederum mit Gl. A5. Mit Gl. A9 ergeben sich die Parameter dk zu:
⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟
⎠
⎞
⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜
⎝
⎛
+ε⋅+ε⋅−εεε+ε⋅ε
εε+ε⋅εεε+ε⋅εεε−ε−ε
=
−
−
−
−
ε−−ε−
−−−−
2N
21N
2N
21N e
1N2Ne
2N
2N2N3N3N
3322
2211
1212
320...0000)(2...0000
000...)(20000...0)(2000...0
M
15
k
k1kk 3
ccdε⋅−
= + ; für 1 < k < N-2
Der noch fehlende Parameter dN-1 ergibt sich mit (A16) zu:
dN-1 = dN-2 Mit Gl. A7 ergeben sich die Parameter bk zu:
k2kkk
k
kk dcb ⋅ε−⋅ε−
εδ
= ; für 1 < k < N-1
Abschließend ist zu bemerken, dass sich bei N Stützstellen die beiden Vektoren Cr
, Yr
aus jeweils N-1 Elementen zusammensetzen und M eine quadratische )1N()1N( −×− - Matrix darstellt. Gegenüber der Randbedingung für natürliche Splines verringert sich somit die Differenzenmatrix M um den Wert 1. (Von den insgesamt N Stützstellen werden offensichtlich nur N-1 zur Bestimmung der Polynomkoeffizienten benötigt. Eine Stützstelle ist überflüssig, was man mit der Bezeichnung „not-a-knot“ zum Ausdruck bringt.) (iii) Splines mit der Randbedingung d3 Mit den vorgegebenen Werten α und β für die Steigung der Polynome in den beiden Randpunkten des Interpolationsintervalls und Gl. A3 ergeben sich mit der Randbedingung (d3) die Beziehungen: α== 11
'1 b)x(P
und
β=−⋅⋅+−⋅⋅+= −−−−−−2
1NN1N1NN1N1NN'
1N )xx(d3)xx(c2b)x(P bzw.:
2N1NN1N1N d3c2b ε⋅⋅+ε⋅⋅+=β −−−
Führt man den Pseudoparameter:
bN = β (A23) ein, so läßt sich der Gültigkeitsbereich von Gl. A7 erweitern zu: 2
kkkkk1k d3c2bb ε⋅⋅+ε⋅⋅+=+ ; für 1 < k < N-1 Hieraus folgt:
( )kkk1k2k
k c2bb3
1d ε⋅⋅−−⋅ε⋅
= + ; für 1 < k < N-1 (A24)
Setzt man Gl. A24 in Gl. A7 ein, so ergibt sich:
16
( )k
kkkk1k3
1kkk c2bbcb
εδ
=ε⋅⋅−−⋅+ε⋅+ + ; für 1 < k < N-1
Stellt man diese Gleichung nach dem Parameter ck um, so erhält man:
⎟⎟⎠
⎞⎜⎜⎝
⎛⋅−−
εδ⋅⋅
ε= + k1k
k
k
kk b2b31c ; für 1 < k < N-1 (A25)
Setzt man weiterhin Gl. A24 in Gl. A9 ein, so ergibt sich:
( )kkk1kk
k1k c2bb1cc ε⋅⋅−−⋅ε
+= ++ ; für 1 < k < N-2
und nach Umformung:
( ) kk1kk
1k cbb1c −−⋅ε
= ++ ; für 1 < k < N-2
Setzt man hierin Gl. A25 ein, so ergibt sich:
( ) ⎟⎟⎠
⎞⎜⎜⎝
⎛⋅−−
εδ⋅⋅
ε−−⋅
ε= +++ k1k
k
k
kk1k
k1k b2b31bb1c
oder
2k
kk
k1k
k1k 3b1b2c
εδ⋅−⋅
ε+⋅
ε= ++ ; für 1 < k < N-2
Substituiert man in Gl. A25 für k → k+1, so erhält man:
⎟⎟⎠
⎞⎜⎜⎝
⎛⋅−−
εδ⋅⋅
ε= ++
+
+
++ 1k2k
1k
1k
1k1k b2b31c ; für 0 < k < N-2
Damit erhält man für 1 < k < N-2:
2k
kk
k1k
k1k2k
1k
1k
1k
3b1b2b2b31εδ⋅−⋅
ε+⋅
ε=⎟⎟
⎠
⎞⎜⎜⎝
⎛⋅−−
εδ⋅⋅
ε ++++
+
+
Stellt man diese Beziehung um, so ergibt sich:
⎟⎟⎠
⎞⎜⎜⎝
⎛εδ
+εδ
⋅=⋅ε
−⋅⎟⎟⎠
⎞⎜⎜⎝
⎛ε
+ε
⋅+⋅ε +
++
++
+2
1k
1k2k
k2k
1k1k
1kkk
k
3b1b112b1; für 1 < k < N-2
Insgesamt ergeben sich die N Gleichungen:
b1 = α
⎟⎟⎠
⎞⎜⎜⎝
⎛εδ
+εδ
⋅=⋅ε
−⋅⎟⎟⎠
⎞⎜⎜⎝
⎛ε
+ε
⋅+⋅ε +
++
++
+2
1k
1k2k
k2k
1k1k
1kkk
k
3b1b112b1; für 1 < k < N-2
17
bN = β Führt man die folgenden Definitionen ein: (i) Parametervektor:
⎟⎟⎟⎟⎟⎟⎟⎟⎟
⎠
⎞
⎜⎜⎜⎜⎜⎜⎜⎜⎜
⎝
⎛
⋅⋅=
−
N
1N
3
2
1
bb
bbb
Br
(ii) Differenzenvektor:
⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟
⎠
⎞
⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜
⎝
⎛
β⋅εδ
+εδ
⋅⋅εδ
+εδ
εδ
+εδ
α⋅
⋅=
−
−
−
−
31
21N
1N2
2N
2N
23
322
2
22
221
1
31
3Yr
(iii) Differenzenmatrix:
so ergibt sich die Matrizengleichung:
YBrr
=⋅M Hieraus ergeben sich die N Polynomparameter bk zu: YB 1
rr⋅= −M
( )( )
( )⎟⎟⎟⎟⎟⎟⎟⎟⎟
⎠
⎞
⎜⎜⎜⎜⎜⎜⎜⎜⎜
⎝
⎛
−+⋅
−+⋅−+⋅
=
−−−− εεεε
εεεε
εεεε
100...00002...0000
...
...000...20000...02000...0001
M
1N1N2N2N
3322
2211
1111
1111
1111
18
Mit den Parametern bk; für 1 < k < N berechen sich die Parameter ck gemäß (A25). Zur Berechnung der Parameter dk wird in Gl. A24 der Parameter ck durch Gl. A25 eliminiert. Damit erhält man:
⎟⎟⎠
⎞⎜⎜⎝
⎛⎟⎟⎠
⎞⎜⎜⎝
⎛⋅−−
εδ⋅⋅−−⋅
ε⋅= ++ k1k
k
kk1k2
kk b2b32bb
31d
woraus sich nach Zusammenfassung für:
⎟⎟⎠
⎞⎜⎜⎝
⎛εδ⋅−+⋅
ε= +
k
kk1k2
kk 2bb1d
ergibt. (iv) Splines mit der Randbedingung d4 Mit Gl. A4 erhält man: 11
'1 b)x(P =
und
21NN1N1NN1N1NN
'1N )xx(d3)xx(c2b)x(P −−−−−− −⋅⋅+−⋅⋅+=
oder
2
1N1N1N1N1NN'
1N d3c2b)x(P −−−−−− ε⋅⋅+ε⋅⋅+= Damit ergibt sich mit der Randbedingung d4, also: )x(P)x(P N
'1N1
'1 −= :
2
1N1N1N1N1N1 d3c2bb −−−−− ε⋅⋅+ε⋅⋅+= (A26) Weiterhin ergibt sich mit A5:
11''
1 c2)x(P ⋅= und
)xx(d6c2)x(P 1NN1N1NN''
1N −−−− −⋅⋅+⋅= bzw. 1N1N1NN
''1N d6c2)x(P −−−− ε⋅⋅+⋅=
Damit ergibt sich mit dem zweiten Teil der Randbedingung d4, also: )x(P)x(P N
''1N1
''1 −= die
Beziehung:
1N1N1N1 d3cc −−− ε⋅⋅+= Umgestellt nach dN-1 ergibt:
19
1N
1N11N 3
ccd−
−− ε⋅
−= (A27)
Aus Gl. A9 folgt:
k
k1kk 3
ccdε⋅−
= + ; für 1 < k < N-2
Faßt man diese Beziehung mit (A27) zusammen, so ergibt sich:
k
k1kk 3
ccdε⋅−
= + ; für 1 < k < N-1 (A28)
mit
cN = c1 (A29) Setzt man diese Beziehung in Gl. A7 ein, so ergibt sich:
k
kk1kk3
1kkk )cc(cb
εδ
=−⋅ε⋅+ε⋅+ + ; für 1 < k < N-1
oder
k
k1kkk3
1k )cc2(b
εδ
=+⋅⋅ε⋅+ + ; für 1 < k < N-1 (A30)
Setzt man (A28) in (A8) ein, so ergibt sich:
)cc(c2bb k1kkkkk1k −⋅ε+ε⋅⋅+= ++ ; für 1 < k < N-2 oder
)cc(bb 1kkkk1k ++ +⋅ε+= ; für 1 < k < N-2 (A31)
Stellt man (A30) nach bk um, so erhält man:
)cc2(b 1kkk31
k
kk ++⋅⋅ε⋅−
εδ
= ; für 1 < k < N-1 (A32)
Substituiert man k → k+1, so ergibt sich:
)cc2(b 2k1k1k31
1k
1k1k +++
+
++ +⋅⋅ε⋅−
εδ
= ; für 0 < k < N-2 (A33)
Setzt man (A32) und (A33) in (A31) ein, so ergibt sich die Beziehung:
)1kkk1kkk31
k
k2k1k1k3
1
1k
1k cc()cc2()cc2( ++++++
+ +⋅ε++⋅⋅ε⋅−εδ
=+⋅⋅ε⋅−εδ
; für 1 < k < N-2
Faßt man zugehörige Terme zusammen, so ergibt sich:
20
⎟⎟⎠
⎞⎜⎜⎝
⎛εδ
−εδ
⋅=⋅ε+⋅ε+ε⋅+⋅ε+
+++++
k
k
1k
1k2k1k1k1kkkk 3cc)(2c ; für 1 < k < N-2 (A34)
Setzt des weiteren (A28) in (A26) ein, so ergibt sich:
)cc(bb N1N1N1N1 +⋅ε+= −−− (A35) Setzt man in (A32) für k = 1, so ergibt sich:
)cc2(b 21131
1
11 +⋅⋅ε⋅−
εδ
=
und setzt man für k = N-1, so ergibt sich:
)cc2(b N1N1N31
1N
1N1N +⋅⋅ε⋅−
εδ
= −−−
−−
Setzt man beide Beziehungen in (A35) ein, so wird:
)cc()cc2()cc2( N1N1NN1N1N31
1N
1N2113
1
1
1 +⋅ε++⋅⋅ε⋅−εδ
=+⋅⋅ε⋅−εδ
−−−−−
−
Faßt man zugehörige Terme zusammen, so ergibt sich:
⎟⎟⎠
⎞⎜⎜⎝
⎛εδ
−εδ
⋅=⋅ε+⋅ε⋅+⋅ε−⋅ε⋅−−
−−−−
1
1
1N
1NN1N1N1N2111 3cc2cc2 (A36)
Insgesamt ergibt sich mit (A34, A36, A29) das folgende Gleichungssystem:
⎟⎟⎠
⎞⎜⎜⎝
⎛εδ
−εδ
⋅=⋅ε+⋅ε+ε⋅+⋅ε+
+++++
k
k
1k
1k2k1k1k1kkkk 3cc)(2c ; für 1 < k < N-2
⎟⎟⎠
⎞⎜⎜⎝
⎛εδ
−εδ
⋅=⋅ε+⋅ε⋅+⋅ε−⋅ε⋅−−
−−−−
1
1
1N
1NN1N1N1N2111 3cc2cc2
0cc N1 =−
21
Mit den Definitionen: (i) Parametervektor:
⎟⎟⎟⎟⎟⎟⎟⎟⎟
⎠
⎞
⎜⎜⎜⎜⎜⎜⎜⎜⎜
⎝
⎛
⋅⋅
=
−
−
N
1N
2N
2
1
ccc
cc
Cr
(ii) Differenzenvektor:
⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟
⎠
⎞
⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜
⎝
⎛
εδ
−εδ
εδ
−εδ
⋅⋅εδ
−εδ
εδ
−εδ
⋅=
−
−
−
−
−
−
0
3Y
1
1
1N
1N
3N
3N
2N
2N
1
1
2
2
1
1
2
2
r
(iii) Differenzenmatrix: läßt sich das Gleichungssystem als Matrizengleichung formulieren:
YCrr
=⋅M Bildet man die inverse Differenzenmatrix M-1, so ergibt sich damit der Parametervektor gemäß:
YC 1rr
⋅= −M Mit den Parametern ck; für 1 < k < N lassen sich dann die Parameter bk gemäß (A32) und dk gemäß (A28) berechnen.
⎟⎟⎟⎟⎟⎟⎟⎟⎟
⎠
⎞
⎜⎜⎜⎜⎜⎜⎜⎜⎜
⎝
⎛
−εε⋅ε−ε⋅−εε+ε⋅ε
εε+ε⋅εεε+ε⋅ε
=
−−
−−−−
100...000120...002
)(2...0000......
000...)(20000...0)(2
M
1N1N11
1N1N2N2N
3322
2211
22
Wie schon erwähnt: Aufgrund der gewählten Randbedingung können die einzelnen Spline-Funktionen, insbesondere in den Randbereichen, erheblich voneinander abweichen wie folgendes Beispiel zeigt.
-1 -0.5 0 0.5 1 1.5 2 2.5 3 3.5 4-4
-3
-2
-1
0
1
2
3
4
↑
y(x)
x →
Stützwerted1d2d3d4
Spline-Interpolation mit verschiedenen Randbedingungen