IDL 中的 SPLINE 函数 允许对数据进行三次插值（至少有 3 个数据点）。虽然 Python 中的 Scipy 库可以使用 UnivariateSpline 和 splrep 函数执行类似的计算，但如果插值设置为三次插值并且我们只有 3 个数据点，这些计算就会中断（SPLINE 不会发生这种情况） ).
这是 SPLINE 在 IDL 中执行操作的简单示例：
> x = [2., 3., 4.]
> y = (x-3)^2
> t = FINDGEN(20)/10.+2
> z = SPLINE(x, y, t)
> print, z
1.00000 0.810662 0.642087 0.493590 0.364684 0.255081 0.164684 0.0935898 0.0420876 0.0106618 0.00000 0.0106618 0.0420876 0.0935898 0.164684 0.255081 0.364684 0.493590 0.642087 0.810662
但是如果我尝试在 Python 中使用 从 scipy.interpolate 导入 splrep, splev
x = np.array([2., 3., 4.])
y = (x-3)**2
t = np.arange(20)/10.+2
z = splev(t, splrep(x, y, k = 3))
或与 从 scipy.interpolate 导入 UnivariateSpline
x = np.array([2., 3., 4.])
y = (x-3)**2
t = np.arange(20)/10.+2
z = UnivariateSpline(x, y, k = 3)(t)
我总是收到此错误消息：
TypeError: m > k must hold
我的理解是，如果 m ≤ k，则当我们必须拟合 m 个数据点时，不可能有唯一的 k 次多项式解。但随之而来的问题是……IDL 中的 SPLINE 如何执行此计算？我怎样才能用Python重现它？
我可以尝试将多项式降低到 k = 2 （二次插值），就像这样
z = splev(t, splrep(x, y, k = 2))
或
z = UnivariateSpline(x, y, k = 2)(t)
两种情况我都会得到：
> print(z)
[1. 0.81 0.64 0.49 0.36 0.25 0.16 0.09 0.04 0.01 0. 0.01 0.04 0.09 0.16 0.25 0.36 0.49 0.64 0.81]
这肯定与 IDL 中的输出类似，除非我们忽略小数点第二位以下的所有内容。
即使在像 SPLINE 那样 k = m 的情况下，如何执行与 Python 中的 SPLINE 相同的计算？
由于这个问题仍然没有答案，我已经将 SPLINE.pro 的 170 行代码翻译成 Python3，如下：
import numpy as np
def spline(X, Y, T, sigma = 1.0):
n = min(len(X), len(Y))
if n <= 2:
print('X and Y must be arrays of 3 or more elements.')
if sigma != 1.0:
sigma = min(sigma, 0.001)
yp = np.zeros(2*n)
delx1 = X[1]-X[0]
dx1 = (Y[1]-Y[0])/delx1
nm1 = n-1
nmp = n+1
delx2 = X[2]-X[1]
delx12 = X[2]-X[0]
c1 = -(delx12+delx1)/(delx12*delx1)
c2 = delx12/(delx1*delx2)
c3 = -delx1/(delx12*delx2)
slpp1 = c1*Y[0]+c2*Y[1]+c3*Y[2]
deln = X[nm1]-X[nm1-1]
delnm1 = X[nm1-1]-x[nm1-2]
delnn = X[nm1]-X[nm1-2]
c1 = (delnn+deln)/(delnn*deln)
c2 = -delnn/(deln*delnm1)
c3 = deln/(delnn*delnm1)
slppn = c3*Y[nm1-2]+c2*Y[nm1-1]+c1*Y[nm1]
sigmap = sigma*nm1/(X[nm1]-X[0])
dels = sigmap*delx1
exps = np.exp(dels)
sinhs = 0.5*(exps-1/exps)
sinhin = 1/(delx1*sinhs)
diag1 = sinhin*(dels*0.5*(exps+1/exps)-sinhs)
diagin = 1/diag1
yp[0] = diagin*(dx1-slpp1)
spdiag = sinhin*(sinhs-dels)
yp[n] = diagin*spdiag
delx2 = X[1:]-X[:-1]
dx2 = (Y[1:]-Y[:-1])/delx2
dels = sigmap*delx2
exps = np.exp(dels)
sinhs = 0.5*(exps-1/exps)
sinhin = 1/(delx2*sinhs)
diag2 = sinhin*(dels*(0.5*(exps+1/exps))-sinhs)
diag2 = np.concatenate([np.array([0]), diag2[:-1]+diag2[1:]])
dx2nm1 = dx2[nm1-1]
dx2 = np.concatenate([np.array([0]), dx2[1:]-dx2[:-1]])
spdiag = sinhin*(sinhs-dels)
for i in range(1, nm1):
diagin = 1/(diag2[i]-spdiag[i-1]*yp[i+n-1])
yp[i] = diagin*(dx2[i]-spdiag[i-1]*yp[i-1])
yp[i+n] = diagin*spdiag[i]
diagin = 1/(diag1-spdiag[nm1-1]*yp[n+nm1-1])
yp[nm1] = diagin*(slppn-dx2nm1-spdiag[nm1-1]*yp[nm1-1])
for i in range(n-2, -1, -1):
yp[i] = yp[i]-yp[i+n]*yp[i+1]
m = len(T)
subs = np.repeat(nm1, m)
s = X[nm1]-X[0]
sigmap = sigma*nm1/s
j = 0
for i in range(1, nm1+1):
while T[j] < X[i]:
subs[j] = i
j += 1
if j == m:
break
if j == m:
break
subs1 = subs-1
del1 = T-X[subs1]
del2 = X[subs]-T
dels = X[subs]-X[subs1]
exps1 = np.exp(sigmap*del1)
sinhd1 = 0.5*(exps1-1/exps1)
exps = np.exp(sigmap*del2)
sinhd2 = 0.5*(exps-1/exps)
exps = exps1*exps
sinhs = 0.5*(exps-1/exps)
spl = (yp[subs]*sinhd1+yp[subs1]*sinhd2)/sinhs+((Y[subs]-yp[subs])*del1+(Y[subs1]-yp[subs1])*del2)/dels
if m == 1:
return spl[0]
else:
return spl
现在一切都按预期工作，结果与 IDL 中的完全相同。
x = np.array([2., 3., 4.])
y = (x-3)**2
t = np.arange(20)/10.+2
z = spline(x, y, t)
> print(z)
[1. 0.81066204 0.64208752 0.49359014 0.36468451 0.25508134
0.16468451 0.09359014 0.04208752 0.01066204 0. 0.01066204
0.04208752 0.09359014 0.16468451 0.25508134 0.36468451 0.49359014
0.64208752 0.81066204]
我得到了一些接近的东西： np.interp(fs_times, fast_spin[:,1], 纪元) 但感谢您的样条代码。它也完全符合我们的 IDL