如何在python中执行三次样条插值？

Question

我有两个列表来描述函数y(x):

x = [0,1,2,3,4,5]
y = [12,14,22,39,58,77]

我想执行三次样条插值,以便在x的域中给出一些值u,例如

u = 1.25

我能找到你(你).

我在SciPy中发现了这个,但我不知道如何使用它.

Answer 1

简短回答:

from scipy import interpolate

def f(x):
    x_points = [ 0, 1, 2, 3, 4, 5]
    y_points = [12,14,22,39,58,77]

    tck = interpolate.splrep(x_points, y_points)
    return interpolate.splev(x, tck)

print(f(1.25))

答案很长:

scipy将样条插值中涉及的步骤分为两个操作,最有可能是计算效率.首先,使用splrep()计算描述样条曲线的系数.splrep返回包含系数的元组数组.其次,将这些系数传递到splev()以实际评估所需点处的样条.这种方法对单次评估来说无疑是不方便的,但由于最常见的用例是从少数函数评估点开始,然后重复使用样条函数来查找插值,这在实践中通常非常有用.

Answer 2

如果没有安装scipy:

import numpy as np
from math import sqrt

def cubic_interp1d(x0, x, y):
    """
    Interpolate a 1-D function using cubic splines.
      x0 : a float or an 1d-array
      x : (N,) array_like
          A 1-D array of real/complex values.
      y : (N,) array_like
          A 1-D array of real values. The length of y along the
          interpolation axis must be equal to the length of x.

    Implement a trick to generate at first step the cholesky matrice L of
    the tridiagonal matrice A (thus L is a bidiagonal matrice that
    can be solved in two distinct loops).

    additional ref: www.math.uh.edu/~jingqiu/math4364/spline.pdf 
    """
    x = np.asfarray(x)
    y = np.asfarray(y)

    # remove non finite values
    # indexes = np.isfinite(x)
    # x = x[indexes]
    # y = y[indexes]

    # check if sorted
    if np.any(np.diff(x) < 0):
        indexes = np.argsort(x)
        x = x[indexes]
        y = y[indexes]

    size = len(x)

    xdiff = np.diff(x)
    ydiff = np.diff(y)

    # allocate buffer matrices
    Li = np.empty(size)
    Li_1 = np.empty(size-1)
    z = np.empty(size)

    # fill diagonals Li and Li-1 and solve [L][y] = [B]
    Li[0] = sqrt(2*xdiff[0])
    Li_1[0] = 0.0
    B0 = 0.0 # natural boundary
    z[0] = B0 / Li[0]

    for i in range(1, size-1, 1):
        Li_1[i] = xdiff[i-1] / Li[i-1]
        Li[i] = sqrt(2*(xdiff[i-1]+xdiff[i]) - Li_1[i-1] * Li_1[i-1])
        Bi = 6*(ydiff[i]/xdiff[i] - ydiff[i-1]/xdiff[i-1])
        z[i] = (Bi - Li_1[i-1]*z[i-1])/Li[i]

    i = size - 1
    Li_1[i-1] = xdiff[-1] / Li[i-1]
    Li[i] = sqrt(2*xdiff[-1] - Li_1[i-1] * Li_1[i-1])
    Bi = 0.0 # natural boundary
    z[i] = (Bi - Li_1[i-1]*z[i-1])/Li[i]

    # solve [L.T][x] = [y]
    i = size-1
    z[i] = z[i] / Li[i]
    for i in range(size-2, -1, -1):
        z[i] = (z[i] - Li_1[i-1]*z[i+1])/Li[i]

    # find index
    index = x.searchsorted(x0)
    np.clip(index, 1, size-1, index)

    xi1, xi0 = x[index], x[index-1]
    yi1, yi0 = y[index], y[index-1]
    zi1, zi0 = z[index], z[index-1]
    hi1 = xi1 - xi0

    # calculate cubic
    f0 = zi0/(6*hi1)*(xi1-x0)**3 + \
         zi1/(6*hi1)*(x0-xi0)**3 + \
         (yi1/hi1 - zi1*hi1/6)*(x0-xi0) + \
         (yi0/hi1 - zi0*hi1/6)*(xi1-x0)
    return f0

if __name__ == '__main__':
    import matplotlib.pyplot as plt
    x = np.linspace(0, 10, 11)
    y = np.sin(x)
    plt.scatter(x, y)

    x_new = np.linspace(0, 10, 201)
    plt.plot(x_new, cubic_interp1d(x_new, x, y))

    plt.show()

Answer 3

如果您安装了scipy版本> = 0.18.0,则可以使用scipy.interpolate中的CubicSpline函数进行三次样条插值.

您可以通过在python中运行以下命令来检查scipy版本:

#!/usr/bin/env python3
import scipy
scipy.version.version

如果你的scipy版本> = 0.18.0,你可以运行以下三维样条插值的示例代码:

#!/usr/bin/env python3

import numpy as np
from scipy.interpolate import CubicSpline

# calculate 5 natural cubic spline polynomials for 6 points
# (x,y) = (0,12) (1,14) (2,22) (3,39) (4,58) (5,77)
x = np.array([0, 1, 2, 3, 4, 5])
y = np.array([12,14,22,39,58,77])

# calculate natural cubic spline polynomials
cs = CubicSpline(x,y,bc_type='natural')

# show values of interpolation function at x=1.25
print('S(1.25) = ', cs(1.25))

## Aditional - find polynomial coefficients for different x regions

# if you want to print polynomial coefficients in form
# S0(0<=x<=1) = a0 + b0(x-x0) + c0(x-x0)^2 + d0(x-x0)^3
# S1(1< x<=2) = a1 + b1(x-x1) + c1(x-x1)^2 + d1(x-x1)^3
# ...
# S4(4< x<=5) = a4 + b4(x-x4) + c5(x-x4)^2 + d5(x-x4)^3
# x0 = 0; x1 = 1; x4 = 4; (start of x region interval)

# show values of a0, b0, c0, d0, a1, b1, c1, d1 ...
cs.c

# Polynomial coefficients for 0 <= x <= 1
a0 = cs.c.item(3,0)
b0 = cs.c.item(2,0)
c0 = cs.c.item(1,0)
d0 = cs.c.item(0,0)

# Polynomial coefficients for 1 < x <= 2
a1 = cs.c.item(3,1)
b1 = cs.c.item(2,1)
c1 = cs.c.item(1,1)
d1 = cs.c.item(0,1)

# ...

# Polynomial coefficients for 4 < x <= 5
a4 = cs.c.item(3,4)
b4 = cs.c.item(2,4)
c4 = cs.c.item(1,4)
d4 = cs.c.item(0,4)

# Print polynomial equations for different x regions
print('S0(0<=x<=1) = ', a0, ' + ', b0, '(x-0) + ', c0, '(x-0)^2  + ', d0, '(x-0)^3')
print('S1(1< x<=2) = ', a1, ' + ', b1, '(x-1) + ', c1, '(x-1)^2  + ', d1, '(x-1)^3')
print('...')
print('S5(4< x<=5) = ', a4, ' + ', b4, '(x-4) + ', c4, '(x-4)^2  + ', d4, '(x-4)^3')

# So we can calculate S(1.25) by using equation S1(1< x<=2)
print('S(1.25) = ', a1 + b1*0.25 + c1*(0.25**2) + d1*(0.25**3))

# Cubic spline interpolation calculus example
    #  https://www.youtube.com/watch?v=gT7F3TWihvk

Answer 4

在我之前的文章中，我编写了一个基于 Cholesky 开发的代码来求解三次算法生成的矩阵。不幸的是，由于平方根函数的原因，它可能在某些点集（通常是非均匀点集）上表现不佳。本着与以前相同的精神，还有另一个想法，即使用 Thomas 算法 (TDMA)（请参阅https://en.wikipedia.org/wiki/Tridiagonal_matrix_algorithm）在其定义循环期间部分求解三对角矩阵。但使用TDMA的条件是至少要求矩阵对角占优。然而，在我们的例子中，它应该是正确的，因为|bi| > |ai| + |ci|, ai = h[i], bi = 2*(h[i]+h[i+1]),ci = h[i+1]且h[i]无条件为正。（参见https://www.cfd-online.com/Wiki/Tridiagonal_matrix_algorithm_- TDMA (Thomas_algorithm)

我再次参考了jingqiu的文档（参见我之前的帖子，不幸的是链接已损坏，但仍然可以在网络缓存中找到它）。

TDMA求解器的优化版本可以描述如下：

def TDMAsolver(a,b,c,d):
""" This function is licenced under: Attribution-ShareAlike 3.0 Unported (CC BY-SA 3.0)
    https://creativecommons.org/licenses/by-sa/3.0/
    Author raphael valentin
    Date 25 Mar 2022
    ref. https://www.cfd-online.com/Wiki/Tridiagonal_matrix_algorithm_-_TDMA_(Thomas_algorithm)
"""
n = len(d)

w = np.empty(n-1,float)
g = np.empty(n, float)

w[0] = c[0]/b[0]
g[0] = d[0]/b[0]

for i in range(1, n-1):
    m = b[i] - a[i-1]*w[i-1]
    w[i] = c[i] / m
    g[i] = (d[i] - a[i-1]*g[i-1]) / m
g[n-1] = (d[n-1] - a[n-2]*g[n-2]) / (b[n-1] - a[n-2]*w[n-2])

for i in range(n-2, -1, -1):
    g[i] = g[i] - w[i]*g[i+1]

return g

当能够获得ai、bi、ci、的每个个体时di，在这两个单循环内组合自然三次样条插值器函数的定义就变得很容易。

def cubic_interpolate(x0, x, y):
""" Natural cubic spline interpolate function
    This function is licenced under: Attribution-ShareAlike 3.0 Unported (CC BY-SA 3.0)
    https://creativecommons.org/licenses/by-sa/3.0/
    Author raphael valentin
    Date 25 Mar 2022
"""
xdiff = np.diff(x)
dydx = np.diff(y)
dydx /= xdiff

n = size = len(x)

w = np.empty(n-1, float)
z = np.empty(n, float)

w[0] = 0.
z[0] = 0.
for i in range(1, n-1):
    m = xdiff[i-1] * (2 - w[i-1]) + 2 * xdiff[i]
    w[i] = xdiff[i] / m
    z[i] = (6*(dydx[i] - dydx[i-1]) - xdiff[i-1]*z[i-1]) / m
z[-1] = 0.

for i in range(n-2, -1, -1):
    z[i] = z[i] - w[i]*z[i+1]

# find index (it requires x0 is already sorted)
index = x.searchsorted(x0)
np.clip(index, 1, size-1, index)

xi1, xi0 = x[index], x[index-1]
yi1, yi0 = y[index], y[index-1]
zi1, zi0 = z[index], z[index-1]
hi1 = xi1 - xi0

# calculate cubic
f0 = zi0/(6*hi1)*(xi1-x0)**3 + \
    zi1/(6*hi1)*(x0-xi0)**3 + \
    (yi1/hi1 - zi1*hi1/6)*(x0-xi0) + \
    (yi0/hi1 - zi0*hi1/6)*(xi1-x0)
return f0

CubicSpline该函数给出的结果与 中的函数/类相同scipy.interpolate，正如我们在下图中看到的那样。

也可以实现一阶和二阶解析导数，其描述方式如下：

f1p = -zi0/(2*hi1)*(xi1-x0)**2 + zi1/(2*hi1)*(x0-xi0)**2 + (yi1/hi1 - zi1*hi1/6) + (yi0/hi1 - zi0*hi1/6)
f2p = zi0/hi1 * (xi1-x0) + zi1/hi1 * (x0-xi0)

然后，很容易验证 f2p[0] 和 f2p[-1] 等于 0，然后插值器函数会产生自然样条。

关于自然样条的附加参考： https://faculty.ksu.edu.sa/sites/default/files/numerical_analysis_9th.pdf#page=167

使用示例：

import matplotlib.pyplot as plt
import numpy as np
x = [-8,-4.19,-3.54,-3.31,-2.56,-2.31,-1.66,-0.96,-0.22,0.62,1.21,3]
y = [-0.01,0.01,0.03,0.04,0.07,0.09,0.16,0.28,0.45,0.65,0.77,1]
x = np.asfarray(x)
y = np.asfarray(y)

plt.scatter(x, y)
x_new= np.linspace(min(x), max(x), 10000)
y_new = cubic_interpolate(x_new, x, y)
plt.plot(x_new, y_new)

from scipy.interpolate import CubicSpline
f = CubicSpline(x, y, bc_type='natural')
plt.plot(x_new, f(x_new), label='ref')
plt.legend()
plt.show()

总之，这个更新的算法将比以前的代码 (O(n)) 执行插值具有更好的稳定性和更快的速度。与 numba 或 cython 结合使用，速度甚至会非常快。最后，它完全独立于 Scipy。重要的是，请注意，与大多数算法一样，有时对数据进行标准化（例如针对大数值或小数值）以获得最佳结果很有用。同样，在这段代码中，我不检查nan值或有序数据。

不管怎样，这次更新对我来说是一个很好的教训，我希望它能对某人有所帮助。如果您发现任何奇怪的情况请告诉我。