Python中数值积分的一种简单方法

逐步编码示例，用于逼近常微分方程组的解

数值积分是一种用于逼近常微分方程 (ODE) 解的技术。数值积分有多种方法；这些在速度、准确性和复杂性方面各不相同。仅举几例，有欧拉法、龙格-库塔法和梯形法则。值得庆幸的是，作为程序员或工程师，您不需要知道这些方法如何工作的确切细节，即可为 ODE 或 ODE 系统获得良好的解决方案。

在深入研究代码之前，让我们简要介绍一下常微分方程是什么以及数值积分是如何工作的。这些是采用以下形式的微分方程：

在这里，y 是 x 的未知函数（我们试图在某些 x 值处逼近该函数），F 是 x、y 和 y 的导数（我们拥有或被给定的）的函数。 n 表示微分方程的阶数或最高导数。函数 F 通常是由物理学家、数学家、科学家等推导出的众所周知的方程。例如，您可以在给定粒子速度函数的情况下求解粒子的位置。在本例中，时间将等价于上述等式中的 x，而速度将等价于 F。

同样，您可能有一个 ODE 系统，您正在寻求解决（或近似）。方程将采用以下形式：

在这里，我们有一组我们感兴趣的 y 值（m 的数量），每个值都有自己的导数和未知解。

涉及 ODE 的最常见问题之一是初始值问题 (IVP)。在这些问题中，您将得到 y 的初始值和上面给出的形式的 ODE。通常，您将使用一阶 ODE。对于一阶 ODE，数值方法采用 ODE 并将其视为斜率。使用这个斜率和 x 上的一个非常小的增量，可以增加 y 的初始值以逼近 y 的下一个值。重复此过程，直到找到您感兴趣的所有 y 值。如果仔细选择 x 的步长，则近似值可以产生非常准确的结果。这也适用于 ODE 系统。

对于本文，我们将求解一组三 (m) 个一阶 (n) 阶常微分方程。这些方程式是任意的，因此请根据需要替换您自己的方程式。

在这里，刻度线表示我们变量的单一时间导数。我们将解决我们的 ODE 系统作为初始值问题。这意味着我们会将 I、L 和 P 的初始值传递给我们的数值积分器。使用这些值，我们的数值积分器将在时间上采取小步长，并使用导数在每个时间步长求解我们的变量。我将在编码部分详细介绍数值积分器及其功能。让我们跳入代码！

Importing Libraries

在这里，我们只是导入运行此代码所需的库和函数。

• SciPy 库中的solve_ivp 用于求解初值问题和数值积分
• matplotlib 中的 pyplot 用于绘制数值积分的结果（定义为 plt 以便于调用）

# Importing Packages
from scipy.integrate import solve_ivp
import matplotlib.pyplot as plt

创建用户定义的函数

下一步是将我们的 ODE 系统定义为用户定义的 Python 函数，称为模型。这将允许 solve_ivp 在数值积分时使用该函数。您将在下面的代码中看到，我们首先从传递给函数的 Y 向量中提取 I、L 和 P 变量（solve_ivp 在幕后传递 t 和 Y）。然后，我们使用这些值来求解我们在本文前面定义的时间导数。我们为这些衍生品创建一个列表并返回它们。

# Model for solve_ivp
def model(t, Y):
    I = Y[0]
    L = Y[1]
    P = Y[2]
    dIdt = 1/2*P
    dLdt = 2*(L*P)**(1/3)
    dPdt = 100/I
    dYdt = [dIdt, dLdt, dPdt]
    return dYdt

数值积分求解

如本文前面所述，我们将使用solve_ivp 对我们的ODE 系统进行数值积分。 solve_ivp 函数有几个必需的参数：用户定义的函数、模型、时间跨度元组 tspan 和一组初始条件 Y0。时间跨度和初始条件是随机选择的，所以随意摆弄这些，看看你会得到什么结果。我们还可以传递额外的可选参数，这将帮助我们定义我们希望我们的解决方案具有什么数值方法和精度。我们将选择 5(4) 阶的 Runge-Kutta 方法和 1/10¹⁰ 的相对容差（相对精度）。 solve_ivp 还有很多其他选项。在此处查看文档。[0]

# Initial Conditions
Y0 = [0.5, 0.2, 0.1]  # [I0, L0, P0]

# Time Span of Interest
tspan = (0, 5)  # (t0, tf)

# Solving ODE
sol = solve_ivp(model, tspan, Y0, method='RK45', rtol=1e-10)
I_sol, L_sol, P_sol = sol.y
time = sol.t

运行数值积分器后，我们会将结果存储在变量 sol 中。然后，我们可以使用点运算符提取变量 I、L 和 P 的时间历史或近似解。此外，我们可以拉出求解变量的时间步长。

Plotting Results

最后，我们可以使用 Pyplot 可视化 solve_ivp 的输出。使用我们收集的时间和时间历史数据，我们可以在同一个图上绘制每个变量。没有标题、轴标签或图例，任何情节都是不完整的，因此我们也将包括这些内容。

# Plotting Results
plt.plot(time, I_sol)
plt.plot(time, L_sol)
plt.plot(time, P_sol)
plt.title('Numerical Integration Results')
plt.xlabel('Time')
plt.ylabel('Solution Values')
plt.legend(['I', 'L', 'P'])
plt.show()

这部分代码的输出显示了使用数值积分器得到的 I、L 和 P 的近似解。

有了这个，我们已经成功地使用solve_ivp解决了我们原来的常微分方程组。如您所见，该过程相对简单；您只需为数值积分器设置参数。一旦你这样做了，它就会一帆风顺。

这篇关于 Python 中 ODE 数值积分的简单方法的文章到此结束。希望您学到了一些东西并将其应用到您的项目中。如果您有任何问题，请发表评论！如果您还没有，请拍手并关注！谢谢！