反转插值以给出与所需插值函数值相关联的变量

Question

我试图使用scipy的插值函数反转插值函数.假设我创建了一个插值函数,

import scipy.interpolate as interpolate
interpolatedfunction = interpolated.interp1d(xvariable,data,kind='cubic')

当我指定一个时,是否有一些函数可以找到x:

interpolatedfunction(x) == a

换句话说,"我希望我的插值函数等于a; xvariable的值是多少,这样我的函数等于a？"

我很欣赏我可以通过一些数值方案来做到这一点,但有一种更简单的方法吗？如果插值函数在x变量中是多值的,该怎么办？

Answer 1

有专门的方法来寻找三次样条的根.最简单的用法是InterpolatedUnivariateSpline对象的.roots()方法:

spl = InterpolatedUnivariateSpline(x, y)
roots = spl.roots()

这个发现所有的根源,而不只是一个,因为一般的解算器(fsolve,brentq,newton,bisect,等)做的.

x = np.arange(20)
y = np.cos(np.arange(20))
spl = InterpolatedUnivariateSpline(x, y)
print(spl.roots())

输出 array([ 1.56669456, 4.71145244, 7.85321627, 10.99554642, 14.13792756, 17.28271674])

但是,您希望将样条曲线等同于某个任意数字a而不是0.一个选项是重建样条曲线(您不能只a从中减去它):

solutions = InterpolatedUnivariateSpline(x, y - a).roots()

请注意,这些都不适用于返回的函数interp1d; 它没有roots方法.对于该函数,使用类似泛型的方法fsolve是一种选择,但是一次只能从中获取一个根.无论如何,interp1d当有更强大的方法进行相同类型的插值时,为什么要使用三次样条？

非面向对象的方式

a可以直接a从样条系数中减去,而不是在从数据中减去后重建样条.这要求我们采用非面向对象的插值方法.具体来说,sproot接受splrep如下准备的tck元组:

tck = splrep(x, y, k=3, s=0)
tck_mod = (tck[0], tck[1] - a, tck[2])
solutions = sproot(tck_mod)    

我不确定在这里搞砸tck是否值得获得,因为无论如何,大部分计算时间都可能在根查找中.但是有其他选择是件好事.

Answer 2

创建插值函数后interp_fn，您可以通过函数的根找到xwhere的值interp_fn(x) == a

interp_fn2 = lambda x: interp_fn(x) - a

有许多选项可以在scipy.optimize. 例如，要使用初始值从 10 开始的牛顿方法：

from scipy import optimize

optimize.newton(interp_fn2, 10)

实际例子

创建一个插值函数，然后找到根 fn(x) == 5

import numpy as np
from scipy import interpolate, optimize

x = np.arange(10)
y = 1 + 6*np.arange(10) - np.arange(10)**2
y2 = 5*np.ones_like(x)
plt.scatter(x,y)
plt.plot(x,y)
plt.plot(x,y2,'k-')
plt.show()

# create the interpolated function, and then the offset
# function used to find the roots

interp_fn = interpolate.interp1d(x, y, 'quadratic')
interp_fn2 = lambda x: interp_fn(x)-5

# to find the roots, we need to supply a starting value
# because there are more than 1 root in our range, we need 
# to supply multiple starting values.  They should be 
# fairly close to the actual root

root1, root2 = optimize.newton(interp_fn2, 1), optimize.newton(interp_fn2, 5)

root1, root2
# returns:
(0.76393202250021064, 5.2360679774997898)

Answer 3

如果您的数据是单调的，您也可以尝试以下操作：

inversefunction = interpolated.interp1d(data, xvariable, kind='cubic')