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 求解常微分方程组。
求解一阶 ode 的刚性或非刚性系统的初值问题
dy/dt = func(y, t, ...) [or func(t, y, ...)]
其中 y 可以是一个向量。
注意
默认情况下，func 的前两个参数的必要顺序与
scipy.integrate.ode类和函数scipy.integrate.solve_ivp使用的系统定义函数中的参数顺序相反。要使用签名func(t, y, ...)的函数，参数 tfirst 必须设置为True。- 参数:
- funccallable(y, t, …) 或 callable(t, y, …)
计算 t 时 y 的导数。如果签名是
callable(t, y, ...)，则参数 tfirst 必须设置为True。func 不得修改 y 中的数据，因为它是由 ODE 求解器内部使用的数据的视图。- y0数组
y 的初始条件（可以是向量）。
- t数组
求解 y 的时间点序列。初始值点应该是该序列的第一个元素。此序列必须是单调递增或单调递减的；允许重复值。
- args元组，可选
传递给函数的额外参数。
- Dfuncallable(y, t, …) 或 callable(t, y, …)
func 的梯度（雅可比矩阵）。如果签名是
callable(t, y, ...)，则参数 tfirst 必须设置为True。Dfun 不得修改 y 中的数据，因为它是由 ODE 求解器内部使用的数据的视图。- col_deriv布尔值，可选
如果 Dfun 定义列向下（更快）的导数，则为 True，否则 Dfun 应定义跨行的导数。
- full_output布尔值，可选
如果返回可选输出的字典作为第二个输出，则为 True
- printmessg布尔值，可选
是否打印收敛消息
- tfirst布尔值，可选
如果为 True，则 func（以及 Dfun，如果给定）的前两个参数必须是
t, y，而不是默认的y, t。在 1.1.0 版本中添加。
- 返回:
- y数组，形状 (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’
所需的整数工作数组的长度
‘mused’
每个成功时间步的方法指示器的向量：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）方法允许的最大阶数。
示例
受重力和摩擦力作用的摆锤的角度 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()
