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

Question

我有一个函数，我想得到它的积分函数，像这样：

也就是说，我需要在 multiple points 处获取值，而不是在点 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)的系列。除了我们没有保存值。

所以我在想，是否可以为集成函数创建保存点，这样当它通过save point 时，该值将被保存并继续到下一个点。因此，通过对20的单一进程，我们也可以得到(-20,-19,-18,-17 ... 18,19,20)的值

Answer 1

对于使用 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()

如果你没有函数f的解析版本，而只有xs和ys作为数据点，那么你可以使用scipy的interp1d函数在数据点之间进行插值并通过以与以前相同的方式对该插值函数进行处理：

from scipy.interpolate import interp1d

f = interp1d(xs, ys)

Answer 2

一个可以实施您概述的策略，方法是仅在短时间间隔（连续 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 值。如所写，此代码在积分范围的开头输出为 0 的函数，即 -20。当然，可以减去对应于 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) 的右侧。