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

How to perform cubic spline interpolation in python?

提问by user112829

I have two lists to describe the function y(x):

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

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

I would like to perform cubic spline interpolation so that given some value u in the domain of x, e.g.

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

u = 1.25

I can find y(u).

我可以找到你（你）。

I found this in SciPy but I am not sure how to use it.

我在 SciPy 中找到了这个，但我不知道如何使用它。

采纳答案by youngmit

Short answer:

简短的回答：

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))

Long answer:

长答案：

scipy separates the steps involved in spline interpolation into two operations, most likely for computational efficiency.

scipy 将样条插值中涉及的步骤分为两个操作，最有可能是为了提高计算效率。

1. The coefficients describing the spline curve are computed, using splrep(). splrep returns an array of tuples containing the coefficients.

2. These coefficients are passed into splev() to actually evaluate the spline at the desired point x(in this example 1.25). xcan also be an array. Calling f([1.0, 1.25, 1.5])returns the interpolated points at 1, 1.25, and 1,5, respectively.

1. 使用 splrep() 计算描述样条曲线的系数。splrep 返回一个包含系数的元组数组。

2. 这些系数被传递到 splev() 以实际评估所需点处的样条x（在本例中为 1.25）。 x也可以是数组。调用f([1.0, 1.25, 1.5])的返回插点1，1.25和1,5分别。

This approach is admittedly inconvenient for single evaluations, but since the most common use case is to start with a handful of function evaluation points, then to repeatedly use the spline to find interpolated values, it is usually quite useful in practice.

这种方法对于单次求值来说无疑是不方便的，但由于最常见的用例是从少数函数求值点开始，然后重复使用样条来查找插值，因此在实践中通常非常有用。

回答by Christophe Roussy

Minimal python3 code:

最小的python3代码：

from scipy import interpolate

if __name__ == '__main__':
    x = [ 0, 1, 2, 3, 4, 5]
    y = [12,14,22,39,58,77]

    # tck : tuple (t,c,k) a tuple containing the vector of knots,
    # the B-spline coefficients, and the degree of the spline.
    tck = interpolate.splrep(x, y)

    print(interpolate.splev(1.25, tck)) # Prints 15.203125000000002
    print(interpolate.splev(...other_value_here..., tck))

Based on comment of cwhy and answer by youngmit

基于 cwhy 的评论和 Youngmit 的回答

回答by raphael valentin

In case, scipy is not installed:

如果未安装 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()

回答by nexayq

If you have scipy version >= 0.18.0 installed you can use CubicSpline function from scipy.interpolate for cubic spline interpolation.

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

You can check scipy version by running following commands in python:

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

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

If your scipy version is >= 0.18.0 you can run following example code for cubic spline interpolation:

如果您的 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