scipy - 如何集成线性插值函数?

Shu*_*iao 5 python numpy scipy

我有一个函数,它是一组相对较大的数据的插值。我使用线性插值,interp1d所以有很多像这样的非平滑尖点。quadscipy的函数会因为尖点而发出警告。我想知道如何在没有警告的情况下进行集成?

谢谢!


感谢所有的答案。在这里,我总结了一些解决方案,以防其他人遇到同样的问题:

  1. 就像@Stelios 所做的那样,用于points避免警告并获得更准确的结果。
  2. 在实践中,点数通常大于 的默认 limit( limit=50) quad,所以我选择quad(f_interp, a, b, limit=2*p.shape[0], points=p)避免所有这些警告。
  3. 如果ab不是数据集的相同起点或终点x,则p可以通过以下方式选择点p = x[where(x>=a and x<=b)]

Ste*_*ios 5

quad接受一个可选参数,称为points。根据文档:

点:(浮点数,整数序列),可选

有界积分区间中的一系列断点,其中可能出现被积函数的局部困难(例如,奇点、不连续点)。序列不必排序。

在您的情况下,“困难”points正是数据点的 x 坐标。下面是一个例子:

import numpy as np 
from scipy.integrate import quad
np.random.seed(123)

# generate random data set 
x = np.arange(0,10)  
y = np.random.rand(10)

# construct a linear interpolation function of the data set 
f_interp = lambda xx: np.interp(xx, x, y)
Run Code Online (Sandbox Code Playgroud)

这是数据点和 的图f_interp在此处输入图片说明

现在调用quad

quad(f_interp,0,9)
Run Code Online (Sandbox Code Playgroud)

返回一系列警告以及

(4.89770017785734, 1.3762838395159349e-05)

如果您提供points参数,即

quad(f_interp,0,9, points = x)
Run Code Online (Sandbox Code Playgroud)

它不发出警告,结果是

(4.8977001778573435, 5.437539505167948e-14)

这也意味着与之前的调用相比,结果的准确性要高得多。


War*_*ser 5

相反interp1d,您可以使用scipy.interpolate.InterpolatedUnivariateSpline. integral(a, b)该插值器具有计算定积分的方法。

这是一个例子:

import numpy as np
from scipy.interpolate import InterpolatedUnivariateSpline
import matplotlib.pyplot as plt


# Create some test data.
x = np.linspace(0, np.pi, 21)
np.random.seed(12345)
y = np.sin(1.5*x) + np.random.laplace(scale=0.35, size=len(x))**3

# Create the interpolator.  Use k=1 for linear interpolation.
finterp = InterpolatedUnivariateSpline(x, y, k=1)

# Create a finer mesh of points on which to compute the integral.
xx = np.linspace(x[0], x[-1], 5*len(x))

# Use the interpolator to compute the integral from 0 to t for each
# t in xx.
qq = [finterp.integral(0, t) for t in xx]

# Plot stuff
p = plt.plot(x, y, '.', label='data')
plt.plot(x, y, '-', color=p[0].get_color(), label='linear interpolation')
plt.plot(xx, qq, label='integral of linear interpolation')
plt.grid()
plt.legend(framealpha=1, shadow=True)
plt.show()
Run Code Online (Sandbox Code Playgroud)

剧情:

阴谋