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

参数:
funccallable(y, t, …) or callable(t, y, …)

计算 y 在 t 处的导数。如果签名为 callable(t, y, ...),则参数 tfirst 必须设置为 Truefunc 不得修改 y 中的数据,因为它是由 ODE 求解器内部使用的数据视图。

y0array

y 的初始条件(可以是向量)。

tarray

求解 y 的时间点序列。初始值点应是此序列的第一个元素。此序列必须是单调递增或单调递减的;允许重复值。

argstuple, optional

传递给函数的额外参数。

Dfuncallable(y, t, …) or callable(t, y, …)

func 的梯度(雅可比矩阵)。如果签名为 callable(t, y, ...),则参数 tfirst 必须设置为 TrueDfun 不得修改 y 中的数据,因为它是由 ODE 求解器内部使用的数据视图。

col_derivbool, optional

如果 Dfun 定义了按列的导数(更快),则为 True;否则 Dfun 应定义按行的导数。

full_outputbool, optional

如果将可选输出字典作为第二个输出返回,则为 True。

printmessgbool, optional

是否打印收敛消息

tfirstbool, optional

如果为 True,func(如果给出 Dfun,也包括 Dfun)的前两个参数必须是 t, y,而不是默认的 y, t

1.1.0 版新增。

返回值:
yarray, shape (len(t), len(y0))

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

infodictdict, only returned if 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, optional

如果这些参数中的任何一个不是 None 或非负数,则假定雅可比矩阵是带状的。它们给出了此带状矩阵中下部和上部非零对角线的数量。对于带状情况,Dfun 应返回一个矩阵,其行包含非零带(从最低对角线开始)。因此,当 ml >=0mu >=0 时,从 Dfun 返回的矩阵 jac 的形状应为 (ml + mu + 1, len(y0))jac 中的数据必须以 jac[i - j + mu, j] 包含第 i 个方程相对于第 j 个状态变量的导数的方式存储。如果 col_deriv 为 True,则必须返回此 jac 的转置。

rtol, atolfloat, optional

输入参数 rtolatol 决定了求解器执行的误差控制。求解器将根据 max-norm of (e / ewt) <= 1 的形式的不等式来控制 y 中估计局部误差向量 e,其中 ewt 是一个正误差权重向量,计算公式为 ewt = rtol * abs(y) + atol。rtol 和 atol 可以是与 y 长度相同的向量,也可以是标量。默认值为 1.49012e-8。

tcritndarray, optional

需要注意积分的临界点(例如,奇点)向量。

h0浮点数, (0: 求解器决定), 可选

第一步尝试的步长。

hmax浮点数, (0: 求解器决定), 可选

允许的最大绝对步长。

hmin浮点数, (0: 求解器决定), 可选

允许的最小绝对步长。

ixpr布尔值, 可选

在方法切换时是否生成额外的打印输出。

mxstep整数, (0: 求解器决定), 可选

t 中每个积分点允许的最大(内部定义的)步数。

mxhnil整数, (0: 求解器决定), 可选

打印的最大消息数。

mxordn整数, (0: 求解器决定), 可选

非刚性 (Adams) 方法允许的最大阶数。

mxords整数, (0: 求解器决定), 可选

刚性 (BDF) 方法允许的最大阶数。

另请参阅

solve_ivp

求解常微分方程组的初值问题

ode

一个基于 VODE 的更面向对象的积分器

quad

用于求解曲线下面积

示例

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

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

其中 bc 是正常数,撇号 (’) 表示导数。要使用 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 以生成解决方案。要将参数 bc 传递给 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