python插值计算

gsgbgxp

已于 2023-03-14 17:40:09 修改

阅读量8.2k

点赞数 12

分类专栏： python 文章标签： python

于 2021-04-26 14:44:23 首次发布

本文链接：https://blog.csdn.net/gsgbgxp/article/details/116145501

版权

python 专栏收录该内容

29 篇文章 7 订阅

订阅专栏

一维插值

常规用法

【参考文章】
python实现各种插值法(数值分析)
python scipy样条插值函数大全（interpolate里interpld函数）
Python：插值interpolate模块
示例代码

# -*-coding:utf-8 -*-
import numpy as np
from scipy import interpolate
import pylab as pl

x=np.linspace(0,10,11)
#x=[ 0. 1. 2. 3. 4. 5. 6. 7. 8. 9. 10.]
y=np.sin(x)
xnew=np.linspace(0,10,101)
pl.plot(x,y,"ro")

for kind in ["nearest","zero","slinear","quadratic","cubic"]:#插值方式
 #"nearest","zero"为阶梯插值
 #slinear 线性插值
 #"quadratic","cubic" 为2阶、3阶B样条曲线插值
 f=interpolate.interp1d(x,y,kind=kind)
 # ‘slinear', ‘quadratic' and ‘cubic' refer to a spline interpolation of first, second or third order)
 ynew=f(xnew)
 pl.plot(xnew,ynew,label=str(kind))
pl.legend(loc="lower right")
pl.show()

处理外推

当插值超出范围的时候如果不加处理会导致报错，这时有如下几个方案。
【参考文章】
插值超出Python范围
 如何使scipy.interpolate给出超出输入范围的外推结果？

在interpolate下设置bounds_error和fill_value

bounds_error：布尔值，可选
如果为True，则任何时候尝试对x范围之外的值进行插值都会引发ValueError（需要进行插值）。如果为False，则为边界值分配fill_value。fill_value默认情况下返回Nan，如果采用外推则设置fill_value=“extrapolate”。

import numpy as np
from scipy import interpolate

x = np.arange(0,10)
y = np.exp(-x/3.0)
f = interpolate.interp1d(x, y, bounds_error=False,fill_value="extrapolate")

print(f(9),np.exp(-9/3),f(9)-np.exp(-9/3))  
#0.049787068367863944 0.049787068367863944 0.0
print(f(11),np.exp(-11/3),f(11)-np.exp(-11/3)) 
#0.049787068367863944 0.025561533206507402 0.024225535161356542

由于外推部分是线性外推，所以对于非线性函数会有些许误差。
当然fill_value也可以是具体的值，比如0，或者边界值y[-1]（但要注意设置边界值的话可能出现左侧溢出却返回右侧边界值的情况，此时应做相应处理，或者采用下面的interp函数。

采用interp

interp是将左侧和右侧的边界值直接作为超出范围的部分，而且使用格式与interpolate有些不同

import numpy as np
from scipy import interp

x = np.arange(0,10)
y = np.exp(-x/3.0)
interp([9,10],x, y)  #两个插值返回结果相同

自编外推插值函数

参考文章如何使scipy.interpolate给出超出输入范围的外推结果？

单变量插值UnivariateSpline

参考文章 Python：插值interpolate模块
UnivariateSpline可以外插值
调用方式如下：

UnivariateSpline(x,y,w=None,bbox=[None,None],k=3,s=None)

x,y是X-Y坐标数组
w是每个数据点的权重值
k为样条曲线的阶数
s为平滑参数。
s=0，样条曲线强制通过所有数据点
s>0,满足 $\sum(w(y-\operatorname{spline}(x)))^{2} \leq s$

s=0强制通过所有数据点的外插值

from scipy import interpolate
import numpy as np
x1=np.linspace(0,10,20)
y1=np.sin(x1)

sx1=np.linspace(0,12,100)
func1=interpolate.UnivariateSpline(x1,y1,s=0)#强制通过所有点
sy1=func1(sx1)

import matplotlib.pyplot as plt
plt.plot(x1,y1,'o')
plt.plot(sx1,sy1)
plt.show()

在这里插入图片描述
也就插值到(0,12)，范围再大就不行了，毕竟插值的专长不在于预测

s>0：不强制通过所有点

import numpy as np
from scipy import interpolate
x2=np.linspace(0,20,200)
y2=np.sin(x2)+np.random.normal(loc=0,scale=1,size=len(x2))*0.2

sx2=np.linspace(0,22,2000)
func2=interpolate.UnivariateSpline(x2,y2,s=8)
sy2=func2(sx2)

import matplotlib.pyplot as plt
plt.plot(x2,y2,'.')
plt.plot(sx2,sy2)
plt.show()

在这里插入图片描述

多维插值

常见用法

后面几个例子中其实比较麻烦的不是插值算法，而是生成数据，而自己常遇到的情况是数表都已经准备好了，直接插值就可以了。
比如拿到一个二维数表z，横坐标就是x索引，纵坐标就是y索引，进而直接用 interp2d 就可以了。

官方文档：scipy.interpolate.interp2d

from scipy import interpolate
f = interpolate.interp2d(x, y, z, kind='linear', copy=True, bounds_error=False, fill_value=None)  # 首先定义插值函数
z1 = f(x1,y1)   # 如果想插值得到 (x1,y1) 处的新数据，则带入f即可

interp2d()

二维插值仍然参考文章如何使scipy.interpolate给出超出输入范围的外推结果？
使用的函数为interp2d()

scipy.interpolate.interp2d(x, y, z, kind='linear', copy=True, bounds_error=False, fill_value=None)

其中kind为可选项，包括{‘linear’, ‘cubic’, ‘quintic’}，默认是‘linear’。

这里有几个注意事项：

The minimum number of data points required along the interpolation axis is (k+1)**2, with k=1 for linear, k=3 for cubic and k=5 for quintic interpolation.
The interpolator is constructed by bisplrep, with a smoothing factor of 0. If more control over smoothing is needed, bisplrep should be used directly.
The coordinates of the data points to interpolate xnew and ynew have to be sorted by ascending order. interp2d is legacy and is not recommended for use in new code. New code should use RegularGridInterpolator instead.

基本用法

例1

import numpy as np 
def func(x,y): 
    return (x+y)*np.exp(-5*(x**2+y**2)) 
x = np.linspace(-1,1,15)
y = np.linspace(-1,1,20)
xx,yy = np.mgrid[-1:1:20j,-1:1:15j]      # 注意这里维度的先后顺序和x、y相反，j代表数据个数
z = func(xx,yy)

from scipy import interpolate
func_inter = interpolate.interp2d(x,y,z,kind='cubic')

xnew = np.linspace(-1,1,100)
ynew = np.linspace(-1,1,100)
znew = func_inter(xnew,ynew)     #xnew, ynew是一维的，输出znew是二维的
xnew,ynew = np.mgrid[-1:1:100j,-1:1:100j]   #统一变成二维，便于下一步画图

import mpl_toolkits.mplot3d
import matplotlib.pyplot as plt
ax=plt.subplot(111,projection='3d')
ax.plot_surface(xnew,ynew,znew)
ax.scatter(xx,yy,z,c='r',marker='^')
plt.show()

mgrid和meshgrid生成数据的不同

首先注意，mgrid生成的是等差数列，这直接限制了其应用范围，而meshgrid可以是任意数组。
关于两个的用法与区别，可参见 https://blog.csdn.net/gsgbgxp/article/details/129534233

例1
使用 mgrid

import numpy as np
from scipy import interpolate

x = np.linspace(-5.01, 5.01, 30)
y = np.linspace(-5.01, 5.01, 40)
xx, yy = np.mgrid[-5.01:5.01:40j,-5.01:5.01:30j] 
z = np.sin(xx**2+yy**3)
f = interpolate.interp2d(x, y, z, kind='cubic')

import matplotlib.pyplot as plt

xnew = np.arange(-5.01, 5.01, 1e-2)
ynew = np.arange(-5.01, 5.01, 1e-2)
znew = f(xnew, ynew)

plt.plot(x, z[0, :], 'ro-', xnew, znew[0, :], 'b-')
plt.show()

运行结果为
在这里插入图片描述

例2
使用 meshgrid

import numpy as np
from scipy import interpolate

x = np.arange(-5.01, 5.01, 0.25)
y = np.arange(-5.01, 5.01, 0.15)
xx, yy = np.meshgrid(x, y)
z = np.sin(xx**2+yy**3)
f = interpolate.interp2d(x, y, z, kind='cubic')

import matplotlib.pyplot as plt

xnew = np.arange(-5.01, 5.01, 1e-2)
ynew = np.arange(-5.01, 5.01, 1e-2)
znew = f(xnew, ynew)

plt.plot(x, z[0, :], 'ro-', xnew, znew[0, :], 'b-')
plt.show()

运行结果为
在这里插入图片描述
例3
发现上面两个例子运行结果不一样，这和mgrid与meshgrid之间的逻辑差别有关。如果想要得到一样的结果，可以修改例2中的z函数表达式。

import numpy as np
from scipy import interpolate

x = np.arange(-5.01, 5.01, 0.25)
y = np.arange(-5.01, 5.01, 0.15)
xx, yy = np.meshgrid(x, y)
z = np.sin(xx**3+yy**2)
f = interpolate.interp2d(x, y, z, kind='cubic')

import matplotlib.pyplot as plt

xnew = np.arange(-5.01, 5.01, 1e-2)
ynew = np.arange(-5.01, 5.01, 1e-2)
znew = f(xnew, ynew)

plt.plot(x, z[0, :], 'ro-', xnew, znew[0, :], 'b-')
plt.show()

运行结果为
在这里插入图片描述

生成间隔不均匀的数据并插值

mgrid只能生成间隔均匀的数据，因此需要采用meshgrid

# 如果x和y间隔不均匀-meshgrid形式

import numpy as np
from scipy import interpolate

x_start = 0
x_stop  = 5
x_num   = 10
x_base  = 1.5
y_start = 0
y_stop  = 5
y_num   = 10
y_base  = 1.5

x = np.logspace(x_start,x_stop,x_num,base=x_base)    # logspace生成等比数列，底为base
y = np.logspace(y_start,y_stop,y_num,base=y_base)
xx, yy = np.meshgrid(x, y)
z = np.sin(xx**2+yy**3)
f = interpolate.interp2d(x, y, z, kind='cubic')

import matplotlib.pyplot as plt

xnew = np.arange(x_base**x_start,x_base**x_stop, 1e-2)
ynew = np.arange(y_base**y_start,y_base**y_stop, 1e-2)
znew = f(xnew, ynew)

plt.plot(x, z[0, :], 'ro-', xnew, znew[0, :], 'b-')
plt.show()

运行结果为
在这里插入图片描述

二维插值Rbf()

Rbf的优点是，排列可以无序，可以不是等距的网格

step1:随机生成点，并计算函数值

import numpy as np
def func(x,y):
    return (x+y)*np.exp(-5*(x**2+y**2))

x=np.random.uniform(low=-1,high=1,size=100)
y=np.random.uniform(low=-1,high=1,size=100)
z=func(x,y)

step2：插值

from scipy import  interpolate
func=interpolate.Rbf(x,y,z,function='multiquadric')
xnew,ynew=np.mgrid[-1:1:100j,-1:1:100j]
znew=func(xnew,ynew)#输入输出都是二维

step3：画图

import mpl_toolkits.mplot3d
import matplotlib.pyplot as plt
ax=plt.subplot(111,projection='3d')
ax.plot_surface(xnew,ynew,znew)
ax.scatter(x,y,z,c='r',marker='^')
plt.show()