Scipy: 高效生成一系列积分（积分函数）

Question

我有一个函数，想要得到它的积分函数，类似于这样：

那就是说，我需要在多个点获取值，而不是在点 x 上获取单个积分值。

例如：

假设我想要在 (-20,20) 范围内获取值。

def f(x):
    return x**2

x_vals  = np.arange(-20, 21, 1)
y_vals =[integrate.nquad(f, [[0, x_val]]) for x_val in x_vals ]

plt.plot(x_vals, y_vals,'-', color = 'r')

问题

在我提供的示例代码中，对于每个点，从头开始做积分。在我的实际代码中，f(x)相当复杂，而且是多次积分，因此运行时间非常慢(Scipy: speed up integration when doing it for the whole surface?)。

我想知道是否有任何有效的方法可以在给定范围内高效地生成Phi(x)。

我的想法：

点Phi(20)处的积分值是从Phi(19)计算得出的，Phi(19)是从Phi(18)等依此类推。因此，当我们得到Phi(20)时，实际上我们也得到了一系列的(-20,-19,-18,-17 ... 18,19,20)。只是我们没有保存这个值。

我在想，是否有可能为一个集成函数创建“保存点”，这样当它经过一个“保存点”时，值就会被保存并继续到下一个点。因此，通过单一的处理过程可以得到20的值，同时也可以得到(-20,-19,-18,-17 ... 18,19,20)的值。

- cqcn1991

Answer 1

通过仅对短间隔（在相邻x值之间）进行积分，然后对结果进行累加，可以实现您概述的策略。像这样：

import numpy as np
import scipy.integrate as si
def f(x):
    return x**2
x_vals = np.arange(-20, 21, 1)
pieces = [si.quad(f, x_vals[i], x_vals[i+1])[0] for i in range(len(x_vals)-1)]
y_vals = np.cumsum([0] + pieces)

这里的pieces是短时间间隔内的积分，将其相加以产生y值。按照编写方式，此代码输出的函数在积分范围开头处-20为0。当然，可以减去与x = 0对应的y值，以获得与您的图表相同的标准化。

话虽如此，拆分和求和过程是不必要的。当您找到f的不定积分时，实际上是在解决微分方程F' = f。SciPy有一个内置的方法odeint来解决这个问题。只需使用它：

import numpy as np
import scipy.integrate as si
def f(x):
    return x**2
x_vals = np.arange(-20, 21, 1)
y_vals = si.odeint(lambda y,x: f(x), 0, x_vals)

输出结果与第一版本基本相同（存在微小的计算误差），但代码更少。使用 lambda y,x: f(x) 的原因是 odeint 的第一个参数必须是一个接受两个参数的函数，即方程 y' = f(y, x) 的右侧。

- user3717023

---

那么这个方法不能应用于由核密度估计生成的函数吗？ - ZK Zhao

我不认为有什么问题，f的定义方式在这里并不重要。 - user3717023

另外，我注意到这是一元情况。我的情况是双变量的（https://dev59.com/5ZPfa4cB1Zd3GeqPFJZh?noredirect=1&lq=1）。我查看了http://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.odeint.html，在这种情况下似乎不适用于我？ - ZK Zhao

Answer 2

对于使用scipy的solve_ivp等效版本的user3717023的答案，您需要记住函数f中x和y的不同排序（与odeint版本不同）。

此外，请记住，您只能计算到一个常数解。因此，您可能需要根据某些给定条件移动结果。在这里的示例中（使用OP提供的函数f（x）= x ^ 2），我将数值解移动，使其通过原点，匹配最简单的解析解F（x）= x ^ 3/3。

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp

def f(x):
    return x**2

xs = np.linspace(-20, 20, 1001)

# This is the integration step:
sol = solve_ivp(lambda x, y: f(x), t_span=(xs[0], xs[-1]), y0=[0], t_eval=xs)

plt.plot(sol.t, sol.t**3/3,                         ls='-',  c='C0', label="analytic: $F(x)=x^3/3$")
plt.plot(sol.t, sol.y[0],                           ls='--', c='C1', label="numeric solution")
plt.plot(sol.t, sol.y[0] - sol.y[0][sol.t.size//2], ls='-.', c='C3', label="shifted solution going through origin")
plt.legend()

如果您没有函数的分析版本，而只有和作为数据点，则可以使用scipy的函数在数据点之间进行插值，并以与之前相同的方式传递该插值函数：

from scipy.interpolate import interp1d

f = interp1d(xs, ys)

- Zaus