scipy.integrate.

odeint#

scipy.integrate.odeint(func, y0, t, args=(), Dfun=None, col_deriv=0, full_output=0, ml=None, mu=None, rtol=None, atol=None, tcrit=None, h0=0.0, hmax=0.0, hmin=0.0, ixpr=0, mxstep=0, mxhnil=0, mxordn=12, mxords=5, printmessg=0, tfirst=False)[源代码][源代码]#

集成一个常微分方程系统。

备注

对于新代码，使用 scipy.integrate.solve_ivp 来求解微分方程。

使用 FORTRAN 库 odepack 中的 lsoda 求解常微分方程组。

解决刚性或非刚性一阶常微分方程系统的初值问题:

dy/dt = func(y, t, ...)  [or func(t, y, ...)]

其中 y 可以是一个向量。

备注

默认情况下，func 的前两个参数的顺序与 scipy.integrate.ode 类和 scipy.integrate.solve_ivp 函数使用的系统定义函数的参数顺序相反。要使用签名 func(t, y, ...) 的函数，必须将参数 tfirst 设置为 True。

参数:

函数callable(y, t, …) 或 callable(t, y, …) 可调用对象: 计算 y 在 t 处的导数。如果签名是 callable(t, y, ...)，那么参数 tfirst 必须设置为 True。func 不得修改 y 中的数据，因为它是 ODE 求解器内部使用的数据的视图。
y0数组: y 的初始条件（可以是一个向量）。
t数组: 要解 y 的时间点序列。初始值点应是此序列的第一个元素。此序列必须单调递增或单调递减；允许重复值。
参数tuple, 可选: 传递给函数的额外参数。
Dfuncallable(y, t, …) 或 callable(t, y, …) 可调用对象: func 的梯度（雅可比矩阵）。如果签名是 callable(t, y, ...)，那么参数 tfirst 必须设置为 True。Dfun 不得修改 y 中的数据，因为它是 ODE 求解器内部使用的数据的视图。
col_derivbool, 可选: 如果 Dfun 按列定义导数（更快），则为 True，否则 Dfun 应按行定义导数。
完整输出bool, 可选: 如果为真，则返回一个可选输出的字典作为第二个输出
printmessgbool, 可选: 是否打印收敛信息
tfirstbool, 可选: 如果为 True，func`（以及如果提供的 `Dfun）的前两个参数必须是 t, y 而不是默认的 y, t。

Added in version 1.1.0.

返回:

yarray, 形状 (len(t), len(y0))

包含在每个所需时间 t 的 y 值的数组，初始值 y0 在第一行。

infodict字典，仅在 full_output == True 时返回

包含额外输出信息的字典

键	意义
‘hu’	每个时间步成功使用的步长向量
‘tcur’	每个时间步长达到的 t 值的向量（将始终至少与输入时间一样大）
`tolsf`	容差比例因子的向量，大于1.0，当检测到请求的精度太高时计算。
‘tsw’	上次方法切换时 t 的值（每个时间步长给出）
‘nst’	累计时间步数
‘nfe’	每个时间步的累积函数评估次数
‘nje’	每个时间步的雅可比评估累积次数
‘nqu’	每个成功步骤的方法顺序向量
‘imxer’	加权局部误差向量（e / ewt）中最大幅值分量的索引，否则在错误返回时为-1
‘lenrw’	所需的双工作数组的长度
‘leniw’	所需的整数工作数组的长度
‘沉思’	每个成功时间步的方法指示器向量：1: adams（非刚性），2: bdf（刚性）

其他参数:

ml, muint, 可选: 如果这些值中任何一个不是 None 或非负数，则假定雅可比矩阵是带状的。这些值给出了带状矩阵中下部和上部非零对角线的数量。对于带状情况，Dfun 应返回一个矩阵，其行包含非零带（从最低对角线开始）。因此，当 ml >=0 或 mu >=0 时，Dfun 返回的矩阵 jac 的形状应为 (ml + mu + 1, len(y0))。jac 中的数据必须存储为 jac[i - j + mu, j] 保存第 i 个方程相对于第 j 个状态变量的导数。如果 col_deriv 为 True，则必须返回此 jac 的转置。
rtol, atolfloat, 可选: 输入参数 rtol 和 atol 决定了求解器执行的误差控制。求解器将根据 max-norm of (e / ewt) <= 1 形式的非等式控制 y 的估计局部误差向量 e，其中 ewt 是正误差权重向量，计算为 ewt = rtol * abs(y) + atol。rtol 和 atol 可以是与 y 长度相同的向量，也可以是标量。默认值为 1.49012e-8。
tcritndarray，可选: 关键点的向量（例如，奇点），在这些点上应特别注意积分。
h0float, (0: 求解器确定), 可选: 在第一步尝试的步长。
hmaxfloat, (0: 求解器确定), 可选: 允许的最大绝对步长。
hminfloat, (0: 求解器确定), 可选: 允许的最小绝对步长。
ixprbool, 可选: 是否在方法切换时生成额外的打印输出。
mxstepint, (0: 由求解器决定), 可选: 每个积分点在 t 中允许的最大（内部定义的）步数
mxhnilint, (0: 由求解器决定), 可选: 打印消息的最大数量。
mxordnint, (0: 由求解器决定), 可选: 非刚性（Adams）方法允许的最大阶数。
mxordsint, (0: 由求解器决定), 可选: 刚性（BDF）方法允许的最大阶数。

参见

solve_ivp: 求解常微分方程组初值问题
ode: 一个基于 VODE 的更面向对象的积分器
quad: 用于求曲线下的面积

示例

受重力和摩擦力作用的摆的角度 theta 的二阶微分方程可以写成:

theta''(t) + b*theta'(t) + c*sin(theta(t)) = 0

其中 b 和 c 是正的常数，而撇号（’）表示导数。要使用 odeint 求解这个方程，我们必须首先将其转换为一组一阶方程。通过定义角速度 omega(t) = theta'(t)，我们得到系统:

theta'(t) = omega(t)
omega'(t) = -b*omega(t) - c*sin(theta(t))

设 y 为向量 [theta, omega]。我们在 Python 中实现这个系统为：

>>> import numpy as np
>>> def pend(y, t, b, c):
...     theta, omega = y
...     dydt = [omega, -b*omega - c*np.sin(theta)]
...     return dydt
...

我们假设常数为 b = 0.25 和 c = 5.0：

>>> b = 0.25
>>> c = 5.0

对于初始条件，我们假设摆锤接近垂直，theta(0) = pi - 0.1，并且最初处于静止状态，因此 omega(0) = 0。那么初始条件向量为

>>> y0 = [np.pi - 0.1, 0.0]

我们将在区间 0 <= t <= 10 内生成一个在 101 个均匀间隔样本的解。因此，我们的时间数组是：

>>> t = np.linspace(0, 10, 101)

调用 odeint 生成解。为了将参数 b 和 c 传递给 pend，我们使用 args 参数将它们传递给 odeint。

>>> from scipy.integrate import odeint
>>> sol = odeint(pend, y0, t, args=(b, c))

解决方案是一个形状为 (101, 2) 的数组。第一列是 theta(t)，第二列是 omega(t)。以下代码绘制了这两个分量。

>>> import matplotlib.pyplot as plt
>>> plt.plot(t, sol[:, 0], 'b', label='theta(t)')
>>> plt.plot(t, sol[:, 1], 'g', label='omega(t)')
>>> plt.legend(loc='best')
>>> plt.xlabel('t')
>>> plt.grid()
>>> plt.show()

../../_images/scipy-integrate-odeint-1.png