Kubische Splines

1. Aufgabenstellung

Gegeben ist eine Menge von ⟪n+1⟫ Punkten ⟪ P_0(t_0,s_0)⟫, ⟪ P_1(t_1,s_1)⟫, …, ⟪P_n(t_n,s_n)⟫ mit ⟪ t_i ⟫ streng monoton steigend.

    ^
    |        n = 4
 s₄ +                           o
 s₃ +                   o       ·
 s₁ +         o         ·       ·
 s₂ |         ·   o     ·       ·
 s₀ +   o     ·   ·     ·       ·
    |   ·     ·   ·     ·       ·
    |   ·     ·   ·     ·       ·
----+---+-----+---+-----+-------+--->
        t0    t1  t2    t3      t4

Gesucht ist eine zusammenhängende Kurve ⟪ f(t) ⟫ mit ⟪ f(t_1)=s_1 ⟫, ⟪ f(t_2)=s_2 ⟫, …, ⟪ f(t_n)=s_n ⟫, die abschnittsweise aus Polynomen dritten Gerades ⟪ f_i(t) ⟫ besteht. Die Funktion soll zweimal stetig differenzierbar sein (weder ⟪ f'(t) ⟫ noch ⟪ f''(t) ⟫ sollen "springen"). Bei einer Trajektorienplanung sind ⟪ t ⟫ die Zeit und ⟪ s ⟫ die Position auf der Achse.

Die Gesamtfunktion:

⟪(1)⟫ ⟪ f(t) = { ( f_0(t) , if t in [ t_0, t_1] "," ), ( f_1(t) , if t in [ t_1, t_2] "," ), ( , vdots ), ( f_i(t) , if t in [ t_i, t_(i+1)] ","), ( , vdots ), ( f_(n-1)(t), if t in [ t_(n-1), t_n] "." ) :}⟫

Die Gesamtfunktion soll stetig sein:

⟪(2)⟫ ⟪{: ( f_(i-1)(t_i) , = f_i(t_i) , = s_i ), ( f_i(t_(i+1)) , = f_(i+1)(t_(i+1)) , = s_(i+1)) :}⟫

Die Gesamtfunktion soll stetig differenzierbar sein:

⟪(3)⟫ ⟪{: ( f_(i-1)'(t_i) , = f_i'(t_i) ), ( f_i'(t_(i+1)) , = f_(i+1)'(t_(i+1)) ) :}⟫

Die Gesamtfunktion soll zwei mal stetig differenzierbar sein:

⟪(4)⟫ ⟪{: ( f_(i-1)''(t_i) , = f_i''(t_i) , = ddot s_i ), ( f_i''(t_(i+1)) , = f_(i+1)''(t_(i+1)) , = ddot s_(i+1) ) :}⟫

2. Teilfunktionen

Wir betrachten eine Teilfunktion ⟪ f'_i(t) ⟫ mit ⟪ t in [ t_i, t_(i+1)] ⟫:

Bestimmung von fi''(t):

Die Funktion ⟪ f_i(t) ⟫ ist ein kubisches Polynom, daher ist deren erste Ableitung ⟪ f'_i(t) ⟫ quadratisch und die zweite Ableitung ⟪ f''_i(t) ⟫ (bei der Trajektorienberechnung die Beschleunigung) linear und durch ihren Wert an zwei Stellen (⟪ t_i ⟫ und ⟪ t_(i+1) ⟫) eindeutig bestimmt.

Wir benutzen die symmetrische Schreibweise:

⟪(5)⟫ ⟪ f''_i(t) = (t_(i+1)-t) · ddot s_i/d_i + (t-t_i) · ddot s_(i+1)/d_i " mit " d_i := t_(i+1) - t_i ⟫

Dabei ist ⟪ ddot s_i = f''_i(t_i) = f''_(i-1)(t_i) ⟫ die Beschleunigung am linken Rand unseres Bereiches und ⟪ ddot s_(i+1) = f''_i(t_(i+1)) = f''_(i+1)(t_(i+1)) ⟫ die Beschleunigung am rechten Rand.

Die Werte der ⟪ ddot s_i ⟫ werden weiter unten bestimmt.

Bestimmung von fi'(t):

Wir nutzen ⟪(5)⟫ und:

⟪int (t_(i+1)-t) dt = t_(i+1)t-1/2 t^2+C =-1/2 (t^2+2t t_(i+1)+t_(i+1)^2)+C+1/2 t_(i+1)^2=-1/2(t_(i+1)-t)^2+C'⟫

⟪int (t-t_i) dt = 1/2 t^2 - t_i t+C = 1/2 (t^2-2t t_i +t_i^2) +C-1/2 t_i^2 = 1/2(t-t_i )^2+C'⟫

und integrieren ⟪ f''_i(t) ⟫:

⟪ {: ( f'_i(t) , = int f''_i(t) dt ), ( , = int ( (t_(i+1)-t) · ddot s_i/d_i + (t-t_i) · ddot s_(i+1)/d_i ) dt ), ( , = int (t_(i+1)-t)dt · ddot s_i/d_i + int (t-t_i)dt · ddot s_(i+1)/d_i ), ( , = (-1/2 (t_(i+1)-t)^2+C_1) · ddot s_i/d_i + (1/2(t-t_i)^2+C_2) · ddot s_(i+1)/d_i ) :} ⟫

⟪(6)⟫ ⟪ {: ( f'_i(t) , = -1/2 (t_(i+1)-t)^2 · ddot s_i/d_i + 1/2(t-t_i)^2 · ddot s_(i+1)/d_i + C ) :} ⟫

Bestimmung von f(t):

Wir nutzen ⟪(6)⟫ und:

⟪ int -(t_(i+1)-t)^2 dt = +1/3 (t_(i+1)-t)^3 + C ⟫

⟪ int +(t - t_i)^2 dt = +1/3 (t - t_i)^3 + C ⟫

und integrieren ⟪ f'_i(t) ⟫:

⟪ {: ( f_i(t) , = int f'_i(t) dt ), ( , = int -1/2 (t_(i+1)-t)^2 · ddot s_i/d_i + 1/2(t-t_i)^2 · ddot s_(i+1)/d_i + c dt ), ( , = int -(t_(i+1)-t)^2 dt · 1/2 ddot s_i/d_i + int (t-t_i)^2 dt · 1/2 ddot s_(i+1)/d_i + int c dt ), ( , = 1/3 (t_(i+1)-t)^3 · 1/2 ddot s_i/d_i + 1/3 (t-t_i)^3 · 1/2 ddot s_(i+1)/d_i + t·p + C ), ( , = 1/6 (t_(i+1)-t)^3 · ddot s_i/d_i + 1/6 (t-t_i)^3 · ddot s_(i+1)/d_i + t·p + C ), ( , = 1/6 (t_(i+1)-t)^3 · ddot s_i/d_i + 1/6 (t-t_i)^3 · ddot s_(i+1)/d_i + (t·p - t_i·p) + (t_i·p + C) ) :} ⟫

⟪(7)⟫ ⟪ {: ( f_i(t) , = 1/6 (t_(i+1)-t)^3 · ddot s_i/d_i + 1/6 (t-t_i)^3 ddot s_(i+1)/d_i + (t-t_i)·p + q ) :} ⟫

Die Konstanten ⟪ ddot s_i ⟫, ⟪ ddot s_(i+1) ⟫, ⟪ p ⟫ und ⟪ q ⟫ sind noch unbestimmt.

Bestimmung von ⟪ p ⟫ und ⟪ q ⟫ aus den Randbedingungen von fi(t):

Zur Erfüllung der Randbedingungen ⟪ f(t_i) = s_i⟫ und ⟪ f(t_(i+1)) = s_(i+1)⟫ ⟪(2)⟫ setzen wir:

⟪ p := (s_(i+1)-s_i)/d_i - 1/6 d_i (ddot s_(i+1) - ddot s_i) ⟫

⟪ q := s_i - 1/6 d_i^2 ddot s_i ⟫

(Herleitung durch Auflösung eines Gleichungssystems) und erhalten damit:

⟪(8)⟫ ⟪ f_i(t) = 1/6 (t_(i+1)-t)^3 · ddot s_i/d_i + 1/6 (t-t_i)^3 · ddot s_(i+1)/d_i + (t-t_i) ( (s_(i+1)-s_i)/d_i - 1/6 d_i(ddot s_(i+1)-ddot s_i)) + s_i - 1/6 d_i^2·ddot s_i ⟫

Zur Überprüfung setzen wir ein:

⟪ {: ( f_i(t_i), = 1/6 ubrace((t_(i+1)-t_i))_(d_i)""^3 · ddot s_i/d_i + ubrace(1/6 (t_i-t_i)^3 · ddot s_(i+1)/d_i)_0 + ubrace((t_i-t_i) ( (s_(i+1)-s_i)/d_i - 1/6 d_i (ddot s_(i+1)-ddot s_i)))_0 + s_i - 1/6d_i^2·ddot s_i ), ( , = 1/6 d_i^3 · ddot s_i/d_i + s_i - 1/6 d_i^2·ddot s_i = s_i " ✓" ) :} ⟫

⟪ {: ( f_i(t_(i+1)), = ubrace( 1/6 (t_(i+1)-t_(i+1))^3 · ddot s_i/d_i )_0 + 1/6 ubrace( (t_(i+1)-t_i) )_(d_i)""^3 · ddot s_(i+1)/d_i + ubrace( (t_(i+1)-t_i) )_(d_i) ((s_(i+1)-s_i)/d_i - 1/6 d_i(ddot s_(i+1)-ddot s_i)) + s_i - 1/6d_i^2·ddot s_i ), ( , = 1/6 d_i^3 · ddot s_(i+1)/d_i + d_i ( (s_(i+1)-s_i)/d_i - 1/6 d_i (ddot s_(i+1)-ddot s_i)) + s_i - 1/6d_i^2·ddot s_i ), ( , = 1/6 d_i^2 · ddot s_(i+1) + s_(i+1) - s_i - 1/6 d_i^2 (ddot s_(i+1)-ddot s_i) + s_i - 1/6d_i^2·ddot s_i ), ( , = 1/6 d_i^2 · ddot s_(i+1) + s_(i+1) - s_i - 1/6 d_i^2·ddot s_(i+1) + 1/6 d_i^2·ddot s_i + s_i - 1/6d_i^2·ddot s_i = s_(i+1) " ✓" ) :} ⟫

Die Konstanten ⟪ ddot s_i ⟫ und ⟪ ddot s_(i+1) ⟫ sind noch unbestimmt.

3. Verbindung der Teilfunktionen

Gefordert ist in ⟪(3)⟫ die stetige Differenzierbarkeit, also ⟪ f_(i-1)'(t_i) = f_i(t_i)⟫.

Ableitung von ⟪(8)⟫ nach ⟪t⟫ ergibt:

⟪(9)⟫ ⟪ f_i'(t) = -1/2(t_(i+1)-t)^2 · ddot s_i /d_i + 1/2(t - t_i)^2 · ddot s_(i+1)/d_i + (s_(i+1)-s_i)/d_i - 1/6 d_i(ddot s_(i+1)-ddot s_i) | i in [0..n-1] ⟫

⟪(10)⟫ ⟪ f_(i-1)'(t) = -1/2(t_i-t)^2 · ddot s_(i-1) /d_(i-1) + 1/2(t - t_(i-1))^2 · ddot s_i/d_(i-1) + (s_i-s_(i-1))/d_(i-1) - 1/6 d_(i-1)(ddot s_i-ddot s_(i-1)) | i in [1..n] ⟫

Einsetzen von ⟪t_i⟫ in ⟪(9)⟫ ergibt:

⟪ {: (f_i'(t_i),=, -1/2 ubrace((t_(i+1)-t_i))_(d_i)""^2·ddot s_i/d_i + ubrace(1/2(t_i-t_i)^2·ddot s_(i+1)/d_i)_0 + (s_(i+1)-s_i)/d_i - 1/6d_i(ddot s_(i+1)-ddot s_i) ), ( ,=, -1/2 d_i^2·ddot s_i/d_i + (s_(i+1)-s_i)/d_i - 1/6 d_i(ddot s_(i+1)-ddot s_i) ), ( ,=, -1/2 d_i ddot s_i + (s_(i+1)-s_i)/d_i - 1/6 d_i ddot s_(i+1) + 1/6 d_i ddot s_i ) :} ⟫

⟪(11)⟫ ⟪ {: (f_i'(t_i),=, -d_i/3 ddot s_i + (s_(i+1)-s_i)/d_i - d_i/6 ddot s_(i+1) | i in [0..n-1]) :} ⟫

Einsetzen von ⟪t_i⟫ in ⟪(10)⟫ ergibt:

⟪ {: (f_(i-1)'(t_i),=, ubrace(-1/2(t_i-t_i)^2·ddot s_(i-1)/d_(i-1))_0 + 1/2 ubrace((t_i-t_(i-1)))_(d_(t-1))""^2·ddot s_i/d_(i-1) + (s_i-s_(i-1))/d_(i-1) - 1/6 d_(i-1)(ddot s_i-ddot s_(i-1)) ), ( ,=, 1/2 d_(i-1)^2·ddot s_i/d_(i-1) + (s_i-s_(i-1))/d_(i-1) - 1/6 d_(i-1)(ddot s_i-ddot s_(i-1)) ), ( ,=, 1/2 d_(i-1)ddot s_i + (s_i-s_(i-1))/d_(i-1) - 1/6 d_(i-1)ddot s_i + 1/6 d_(i-1)ddot s_(i-1) ) :} ⟫

⟪(12)⟫ ⟪ {: (f_(i-1)'(t_i),=, d_(i-1)/3 ddot s_i + (s_i-s_(i-1))/d_(i-1) + d_(i-1)/6 ddot s_(i-1) | i in [1..n]) :} ⟫

Durch Gleichsetzung von ⟪(12)⟫ und ⟪(11)⟫ erhalten wir:

⟪ {: ( d_(i-1)/3 ddot s_i + (s_i-s_(i-1))/d_(i-1) + d_(i-1)/6 ddot s_(i-1) ,=, -d_i/3 ddot s_i + (s_(i+1)-s_i)/d_i - d_i/6 ddot s_(i+1) ), ( d_(i-1)/3 ddot s_i + (s_i-s_(i-1))/d_(i-1) + d_(i-1)/6 ddot s_(i-1) ,=, -d_i/3 ddot s_i + (s_(i+1)-s_i)/d_i - d_i/6 ddot s_(i+1) ), ( d_(i-1)/6 ddot s_(i-1) + ubrace(d_(i-1)/3 ddot s_i + d_i/3 ddot s_i) + d_i/6 ddot s_(i+1),=, (s_(i+1)-s_i)/d_i - (s_i-s_(i-1))/d_(i-1) ) :} ⟫

⟪(13)⟫ ⟪ d_(i-1)/6 ddot s_(i-1) + (d_(i-1)+d_i)/3 ddot s_i + d_i/6 ddot s_(i+1) = (s_(i+1)-s_i)/d_i - (s_i-s_(i-1))/d_(i-1) | i in [ 1 … n-1 ] ⟫

Wir können diese Gleichungen als Gleichungssytem darstellen:

⟪(14)⟫ ⟪ [ ( ? ,? , , , , , ), ( d_0/6,(d_0+d_1)/3,d_1/6 , , , , ), ( ,ddots ,ddots ,ddots , , , ), ( , ,d_(i-1)/6,(d_(i-1)+d_i)/3,d_i/6 , , ), ( , , ,ddots ,ddots ,ddots , ), ( , , , ,d_(n-2)/6,(d_(n-2)+d_(n-1))/3,d_(n-1)/6 ), ( , , , , ,? ,? ) ] * [ ( ddot s_0 ), ( ddot s_1 ), ( vdots ), ( ddot s_i ), ( vdots ), ( ddot s_(n-1) ), ( ddot s_n ) ] = [ ( ? ), ( (s_2 - s_1) / d_i - (s_1 - s_0) / d_0 ), ( vdots ), ( (s_(i+1) - s_i) / d_i - (s_i - s_(i-1)) / d_(i-1) ), ( vdots ), ( (s_n - s_(n-1)) / d_(n-1) - (s_(n-1) - s_(n-2)) / d_(n-2) ), ( ? ) ] ⟫

4. Randbedingungen

Das Gleichungssystem hat ⟪n+1⟫ Unbekannte, aber nur ⟪n-1⟫ Gleichungen, brauchen also noch zwei zusätzliche Bedingungen.

Dazu können wir Werte von ⟪f_i'(t)⟫ oder ⟪f_i''(t)⟫ an zwei Stellen festlegen.

4.1. Freier Rand (natürliche Randbedingungen)

Wir legen die Krümmung der Randpunkte auf ⟪0⟫ fest (⟪f''(t_0)=0⟫, ⟪f''(t_n)=0⟫) die Randpunkte übertragen also keine "Krümmungsenergie" auf die Funktion. Damit erhalten wir die zusätzlichen Gleichungszeilen:

⟪{: ( 1 * ddot s_0 = 0 ), ( 1 * ddot s_n = 0 ) :}⟫

Diese Zeilen vollständigen das Gleichungssystem zu:

⟪(15)⟫ ⟪ [ ( 1 ,0 , , , , , ), ( d_0/6,(d_0+d_1)/3,d_1/6 , , , , ), ( ,ddots ,ddots ,ddots , , , ), ( , ,d_(i-1)/6,(d_(i-1)+d_i)/3,d_i/6 , , ), ( , , ,ddots ,ddots ,ddots , ), ( , , , ,d_(n-2)/6,(d_(n-2)+d_(n-1))/3,d_(n-1)/6 ), ( , , , , ,0 ,1 ) ] * [ ( ddot s_0 ), ( ddot s_1 ), ( vdots ), ( ddot s_i ), ( vdots ), ( ddot s_(n-1) ), ( ddot s_n ) ] = [ ( 0 ), ( (s_2 - s_1) / d_i - (s_1 - s_0) / d_0 ), ( vdots ), ( (s_(i+1) - s_i) / d_i - (s_i - s_(i-1)) / d_(i-1) ), ( vdots ), ( (s_n - s_(n-1)) / d_(n-1) - (s_(n-1) - s_(n-2)) / d_(n-2) ), ( 0 ) ] ⟫

Die Steigungen ⟪f'(t_0)⟫ und ⟪f'(t_n)⟫ an den Randpunkten ergeben sich, sie sind im allgemeinen ungleich null.

4.2. Eingespannter Rand (hermite Randbedingungen)

Wir setzen die Steigung am Randpunkt ⟪t_0⟫ fest mit ⟪f_0'(t_0)=dot s_0⟫.

⟪ dot s_0 = f_0'(t_0) = -1/2 ubrace((t_1-t_0))_(d_0)""^2 · ddot s_0/d_0 + ubrace(1/2(t_0-t_0)^2 · ddot s_1/d_0)_0 + (s_1-s_0)/d_0 - 1/6 d_0(ddot s_1-ddot s_0) ⟫

⟪ dot s_0 = ubrace(-d_0/2 · ddot s_0) + (s_1-s_0)/d_0 - d_0/6 · ddot s_1 ubrace( +d_0/6 ddot s_0) ⟫

⟪ dot s_0 = -d_0/3 · ddot s_0 + (s_1-s_0)/d_0 - d_0/6 · ddot s_1 ⟫

⟪(16)⟫ ⟪ d_0/3 · ddot s_0 + d_0/6 · ddot s_1 = (s_1-s_0)/d_0 - dot s_0 ⟫

Wir setzen die Steigung am Randpunkt ⟪t_n⟫ fest mit ⟪f_(n-1)'(t_n)=dot s_n⟫.

⟪ dot s_n = f_(n-1)'(t_n) = ubrace(-1/2(t_n-t_n)^2 · ddot s_(n-1)/d_(n-1))_0 + 1/2 ubrace((t_n-t_(n-1)))_(d_(n-1))""^2 · ddot s_n/d_(n-1) + (s_n-s_(n-1))/d_(n-1) - 1/6 d_(n-1)(ddot s_n-ddot s_(n-1)) ⟫

⟪ dot s_n = ubrace(d_(n-1)/2 · ddot s_n) + (s_n-s_(n-1))/d_(n-1) ubrace(-d_(n-1)/6 ddot s_n) + d_(n-1)/6 ddot s_(n-1) ⟫

⟪ dot s_n = d_(n-1)/3 · ddot s_n + (s_n-s_(n-1))/d_(n-1) + d_(n-1)/6 ddot s_(n-1) ⟫

⟪(17)⟫ ⟪ d_(n-1)/6 ddot s_(n-1) + d_(n-1)/3 · ddot s_n = dot s_n - (s_n-s_(n-1))/d_(n-1) ⟫

Mit ⟪(16)⟫ und ⟪(17)⟫ ergibt sich das vollständige Gleichungssystem:

⟪(18)⟫ ⟪ [ ( d_0/3,d_0/6 , , , , , ), ( d_0/6,(d_0+d_1)/3,d_1/6 , , , , ), ( ,ddots ,ddots ,ddots , , , ), ( , ,d_(i-1)/6,(d_(i-1)+d_i)/3,d_i/6 , , ), ( , , ,ddots ,ddots ,ddots , ), ( , , , ,d_(n-2)/6,(d_(n-2)+d_(n-1))/3,d_(n-1)/6 ), ( , , , , ,d_(n-1)/6 ,d_(n-1)/3 ) ] * [ ( ddot s_0 ), ( ddot s_1 ), ( vdots ), ( ddot s_i ), ( vdots ), ( ddot s_(n-1) ), ( ddot s_n ) ] = [ ( (s_1-s_0)/d_0 - dot s_0 ), ( (s_2 - s_1) / d_i - (s_1 - s_0) / d_0 ), ( vdots ), ( (s_(i+1) - s_i) / d_i - (s_i - s_(i-1)) / d_(i-1) ), ( vdots ), ( (s_n - s_(n-1)) / d_(n-1) - (s_(n-1) - s_(n-2)) / d_(n-2) ), ( dot s_n - (s_n-s_(n-1))/d_(n-1)) ] ⟫

5. Algorithmus

⚠ Zusammenfassenden Lösungsalgorithmus zu ergänzen