插值

In [1]:

  1. import numpy as np
  2. import matplotlib.pyplot as plt
  3. %matplotlib inline

设置 Numpy 浮点数显示格式：

In [2]:

  1. np.set_printoptions(precision=2, suppress=True)

从文本中读入数据，数据来自 http://kinetics.nist.gov/janaf/html/C-067.txt ，保存为结构体数组：

In [3]:

  1. data = np.genfromtxt("JANAF_CH4.txt", 
  2.                   delimiter="\t", # TAB 分隔
  3.                   skiprows=1,     # 忽略首行
  4.                   names=True,     # 读入属性
  5.                   missing_values="INFINITE",  # 缺失值
  6.                   filling_values=np.inf)      # 填充缺失值

显示部分数据：

In [4]:

  1. for row in data[:7]:
  2.     print "{}\t{}".format(row['TK'], row['Cp'])
  3. print "...\t..."
  1. 0.0    0.0
  2. 100.0    33.258
  3. 200.0    33.473
  4. 250.0    34.216
  5. 298.15    35.639
  6. 300.0    35.708
  7. 350.0    37.874
  8. ...    ...

绘图：

In [5]:

  1. p = plt.plot(data['TK'], data['Cp'], 'kx')
  2. t = plt.title("JANAF data for Methane $CH_4$")
  3. a = plt.axis([0, 6000, 30, 120])
  4. x = plt.xlabel("Temperature (K)")
  5. y = plt.ylabel(r"$C_p$ ($\frac{kJ}{kg K}$)")

04.02 插值 - 图1

插值

假设我们要对这组数据进行插值。

先导入一维插值函数 interp1d：

  1. interp1d(x, y)

In [6]:

  1. from scipy.interpolate import interp1d

In [7]:

  1. ch4_cp = interp1d(data['TK'], data['Cp'])

interp1d 的返回值可以像函数一样接受输入，并返回插值的结果。

单个输入值，注意返回的是数组：

In [8]:

  1. ch4_cp(382.2)

Out[8]:

  1. array(39.565144000000004)

输入数组，返回的是对应的数组：

In [9]:

  1. ch4_cp([32.2,323.2])

Out[9]:

  1. array([ 10.71,  36.71])

默认情况下，输入值要在插值允许的范围内，否则插值会报错：

In [10]:

  1. ch4_cp(8752)
  1. ---------------------------------------------------------------------------
  2. ValueError                                Traceback (most recent call last)
  3. <ipython-input-10-5d727af9aa33> in <module>()
  4. ----> 1 ch4_cp(8752)
  5.  
  6. d:\Miniconda\lib\site-packages\scipy\interpolate\polyint.pyc in __call__(self, x)
  7.      77         """
  8.      78         x, x_shape = self._prepare_x(x)
  9. ---> 79         y = self._evaluate(x)
  10.      80         return self._finish_y(y, x_shape)
  11.      81 
  12.  
  13. d:\Miniconda\lib\site-packages\scipy\interpolate\interpolate.pyc in _evaluate(self, x_new)
  14.     496         #    The behavior is set by the bounds_error variable.
  15.     497         x_new = asarray(x_new)
  16. --> 498         out_of_bounds = self._check_bounds(x_new)
  17.     499         y_new = self._call(self, x_new)
  18.     500         if len(y_new) > 0:
  19.  
  20. d:\Miniconda\lib\site-packages\scipy\interpolate\interpolate.pyc in _check_bounds(self, x_new)
  21.     526                 "range.")
  22.     527         if self.bounds_error and above_bounds.any():
  23. --> 528             raise ValueError("A value in x_new is above the interpolation "
  24.     529                 "range.")
  25.     530 
  26.  
  27. ValueError: A value in x_new is above the interpolation range.

但我们可以通过参数设置允许超出范围的值存在：

In [11]:

  1. ch4_cp = interp1d(data['TK'], data['Cp'], 
  2.                   bounds_error=False)

不过由于超出范围，所以插值的输出是非法值：

In [12]:

  1. ch4_cp(8752)

Out[12]:

  1. array(nan)

可以使用指定值替代这些非法值：

In [13]:

  1. ch4_cp = interp1d(data['TK'], data['Cp'], 
  2.                   bounds_error=False, fill_value=-999.25)

In [14]:

  1. ch4_cp(8752)

Out[14]:

  1. array(-999.25)

线性插值

interp1d 默认的插值方法是线性，关于线性插值的定义，请参见：

应用线性插值：

In [15]:

  1. T = np.arange(100,355,5)
  2. plt.plot(T, ch4_cp(T), "+k")
  3. p = plt.plot(data['TK'][1:7], data['Cp'][1:7], 'ro', markersize=8)

04.02 插值 - 图2

其中红色的圆点为原来的数据点，黑色的十字点为对应的插值点，可以明显看到，相邻的数据点的插值在一条直线上。

多项式插值

我们可以通过 kind 参数来调节使用的插值方法，来得到不同的结果：

  • nearest 最近邻插值
  • zero 0阶插值
  • linear 线性插值
  • quadratic 二次插值
  • cubic 三次插值
  • 4,5,6,7 更高阶插值 最近邻插值：

In [16]:

  1. cp_ch4 = interp1d(data['TK'], data['Cp'], kind="nearest")
  2. p = plt.plot(T, cp_ch4(T), "k+")
  3. p = plt.plot(data['TK'][1:7], data['Cp'][1:7], 'ro', markersize=8)

04.02 插值 - 图3

0阶插值：

In [17]:

  1. cp_ch4 = interp1d(data['TK'], data['Cp'], kind="zero")
  2. p = plt.plot(T, cp_ch4(T), "k+")
  3. p = plt.plot(data['TK'][1:7], data['Cp'][1:7], 'ro', markersize=8)

04.02 插值 - 图4

二次插值：

In [18]:

  1. cp_ch4 = interp1d(data['TK'], data['Cp'], kind="quadratic")
  2. p = plt.plot(T, cp_ch4(T), "k+")
  3. p = plt.plot(data['TK'][1:7], data['Cp'][1:7], 'ro', markersize=8)

04.02 插值 - 图5

三次插值：

In [19]:

  1. cp_ch4 = interp1d(data['TK'], data['Cp'], kind="cubic")
  2. p = plt.plot(T, cp_ch4(T), "k+")
  3. p = plt.plot(data['TK'][1:7], data['Cp'][1:7], 'ro', markersize=8)

04.02 插值 - 图6

事实上，我们可以使用更高阶的多项式插值，只要将 kind 设为对应的数字即可：

四次多项式插值：

In [20]:

  1. cp_ch4 = interp1d(data['TK'], data['Cp'], kind=4)
  2. p = plt.plot(T, cp_ch4(T), "k+")
  3. p = plt.plot(data['TK'][1:7], data['Cp'][1:7], 'ro', markersize=8)

04.02 插值 - 图7

可以参见：

In [21]:

  1. from scipy.interpolate import interp2d, interpnd

其使用方法与一维类似。

径向基函数

关于径向基函数，可以参阅：

\Phi(x,c) = \Phi(|x-c|)

In [22]:

  1. x = np.linspace(-3,3,100)

常用的径向基（RBF）函数有：

高斯函数：

In [23]:

  1. plt.plot(x, np.exp(-1 * x **2))
  2. t = plt.title("Gaussian")

04.02 插值 - 图8

Multiquadric 函数：

In [24]:

  1. plt.plot(x, np.sqrt(1 + x **2))
  2. t = plt.title("Multiquadric")

04.02 插值 - 图9

Inverse Multiquadric 函数：

In [25]:

  1. plt.plot(x, 1. / np.sqrt(1 + x **2))
  2. t = plt.title("Inverse Multiquadric")

04.02 插值 - 图10

径向基函数插值

对于径向基函数，其插值的公式为：

f(x) = \sum_j n_j \Phi(|x-x_j|)

我们通过数据点 $x_j$ 来计算出 $n_j$ 的值，来计算 $x$ 处的插值结果。

In [26]:

  1. from scipy.interpolate.rbf import Rbf

使用 multiquadric 核的：

In [27]:

  1. cp_rbf = Rbf(data['TK'], data['Cp'], function = "multiquadric")
  2. plt.plot(data['TK'], data['Cp'], 'k+')
  3. p = plt.plot(data['TK'], cp_rbf(data['TK']), 'r-')

04.02 插值 - 图11

使用 gaussian 核：

In [28]:

  1. cp_rbf = Rbf(data['TK'], data['Cp'], function = "gaussian")
  2. plt.plot(data['TK'], data['Cp'], 'k+')
  3. p = plt.plot(data['TK'], cp_rbf(data['TK']), 'r-')

04.02 插值 - 图12

使用 nverse_multiquadric 核：

In [29]:

  1. cp_rbf = Rbf(data['TK'], data['Cp'], function = "inverse_multiquadric")
  2. plt.plot(data['TK'], data['Cp'], 'k+')
  3. p = plt.plot(data['TK'], cp_rbf(data['TK']), 'r-')

04.02 插值 - 图13

不同的 RBF 核的结果也不同。

高维 RBF 插值

In [30]:

  1. from mpl_toolkits.mplot3d import Axes3D

三维数据点：

In [31]:

  1. x, y = np.mgrid[-np.pi/2:np.pi/2:5j, -np.pi/2:np.pi/2:5j]
  2. z = np.cos(np.sqrt(x**2 + y**2))

In [32]:

  1. fig = plt.figure(figsize=(12,6))
  2. ax = fig.gca(projection="3d")
  3. ax.scatter(x,y,z)

Out[32]:

  1. <mpl_toolkits.mplot3d.art3d.Path3DCollection at 0x176b4da0>

04.02 插值 - 图14

3维 RBF 插值：

In [33]:

  1. zz = Rbf(x, y, z)

In [34]:

  1. xx, yy = np.mgrid[-np.pi/2:np.pi/2:50j, -np.pi/2:np.pi/2:50j]
  2. fig = plt.figure(figsize=(12,6))
  3. ax = fig.gca(projection="3d")
  4. ax.plot_surface(xx,yy,zz(xx,yy),rstride=1, cstride=1, cmap=plt.cm.jet)

Out[34]:

  1. <mpl_toolkits.mplot3d.art3d.Poly3DCollection at 0x176e5c50>

04.02 插值 - 图15

原文: https://nbviewer.jupyter.org/github/lijin-THU/notes-python/blob/master/04-scipy/04.02-interpolation-with-scipy.ipynb