ich arbeite seit tagen an einem spline-problem; und ich bekomme es einfach nicht hin!
Kurz: Ich habe eine Punktliste (x,y) welche einen deformierten Kreis mit wenigen
Punkte beschreiben. Ich würde gerne diese Punkte mit einem Spline interpolieren.
Problem: Da die Punkte einen Kreis-beschreiben, soll der Spline an den Randpunkten auch stetig sein.
Ich habe schon diverse Spline funktionen wie interp1d oder splrep probiert, doch am Ende
sehe ich immer eine Unstetigkeitsstelle an den Rändern.
Ich denke ich brauche für meine "geschlossene" interpolation eine periodische Spline-Funktion.
Kenn hier jemand eine solche Funktion?
Ich habe schon probiert mir den Spline selbst zu berechnen - bin bisher aber gescheitert.
Zur basis meiner Berechnungen habe ich diesen Buchauszug benutzt (seite 109): http://www.home.hs-karlsruhe.de/~weth00 ... /kap20.pdf
Doch ich verstehe nicht, wie man genau nach der Lösung des Gleichungssystems auf die Splinekoeffizienten kommt. Kennt sich damit jemand vielleicht aus?
Der vollständigkeitshalber poste hier meinem Versuch für die periodische Spline funktion:
Code: Alles auswählen
def gauss(A): # abkopiert von einem beispiel
n = len(A)
for i in range(0, n):
# Search for maximum in this column
maxEl = abs(A[i][i])
maxRow = i
for k in range(i+1, n):
if abs(A[k][i]) > maxEl:
maxEl = abs(A[k][i])
maxRow = k
# Swap maximum row with current row (column by column)
for k in range(i, n+1):
tmp = A[maxRow][k]
A[maxRow][k] = A[i][k]
A[i][k] = tmp
# Make all rows below this one 0 in current column
for k in range(i+1, n):
c = -A[k][i]/A[i][i]
for j in range(i, n+1):
if i == j:
A[k][j] = 0
else:
A[k][j] += c * A[i][j]
# Solve equation Ax=b for an upper triangular matrix A
x = [0 for i in range(n)]
for i in range(n-1, -1, -1):
x[i] = A[i][n]/A[i][i]
for k in range(i-1, -1, -1):
A[k][n] -= A[k][i] * x[i]
return x
def perispline(x,y):
l = len(x)
num_eq = l-1
em_rows = [0]*(num_eq -1 )
h = [x[i+1] - x[i] for i in range(l-1)]
# first row
r = [0.] * (num_eq -1 + 1)
r[0] = 2*(h[0] + h[1])
r[1] = h[1]
r[len(r) -2] = h[0]
em_rows[0] = r
# body rows
i = 0
while i < len(em_rows)-2:
r = [0.] * (num_eq -1 + 1)
r[i+0] = h[i+1]
r[i+1] = 2*(h[i+1]+h[i+2])
r[i+2] = h[i+2]
em_rows[i+1] = r
i += 1
# last row
r = [0.] * (num_eq -1 + 1)
r[0] = h[0]
r[len(r)-3] = h[len(h) - 2]
r[len(r)-2] = 2 * (h[len(h) - 2] + h[len(h) - 1])
em_rows[len(em_rows)-1] = r
# result column
i = 0
while(i < len(em_rows)):
em_rows[i][len(em_rows[i])-1] = 3. / h[i+1] * (y[i+2] - y[i+1]) - 3. / h[i+0] * (y[i+1] - y[i+0])
i += 1
# solve equation system for c-coeff
gr = gauss(em_rows)
# ich vermute das es bis hier richtig ist....
# coeffs
a = [0.]*num_eq
b = [0.]*num_eq
c = [0.]*num_eq
d = [0.]*num_eq
# set c and a coeffs
i = 0
while i < num_eq:
if i == 0:
c[0] = gr[-1]
else:
c[i] = gr[i-1]
a[i] = y[i-1]
i += 1
i = 0
while i < num_eq:
b[i] = (a[(i+1)%num_eq] - a[i])/h[i] - h[i]/3.*(c[(i+1)%num_eq] + 2*c[i])
d[i] = (c[(i+1)%num_eq] - 2*c[i]) / (3.*h[i])
i += 1
def pericalc(x_value):
#find interval
i = 0
while i < len(x)-1:
if (x_value >= x[i] and x_value <= x[i+1]):
break
i += 1
return a[i] + b[i]*(x_value - x[i]) + c[i]*(x_value - x[i])**2 + d[i]*(x_value - x[i])**3
return pericalc
f = perispline([0.,1.,2.,3.,4.,5.,6.], [1.,2.,1.,2.,1.,2., 1.])
print f(0)
print f(1)
print f(2)
print f(3)
print f(4)
import matplotlib.pyplot as plt
import numpy as np
x = np.linspace(0,4,100)
y = []
for xx in x:
y.append(f(xx))
plt.plot(x,y)
plt.show()
