SciPy-1-12-中文文档-四-

SciPy 1.12 中文文档(四)

原文:docs.scipy.org/doc/scipy-1.12.0/index.html

scipy.fftpack.dstn

原文:docs.scipy.org/doc/scipy-1.12.0/reference/generated/scipy.fftpack.dstn.html#scipy.fftpack.dstn

scipy.fftpack.dstn(x, type=2, shape=None, axes=None, norm=None, overwrite_x=False)

返回沿指定轴进行的多维离散正弦变换。

参数:

x 类似数组。

输入数组。

type {1, 2, 3, 4},可选。

DST 的类型(见注释)。默认类型为 2。

shape 整数或整数数组或 None,可选。

结果的形状。如果 shapeaxes(见下文)均为 None,则 shapex.shape;如果 shape 为 None 但 axes 不为 None,则 shapenumpy.take(x.shape, axes, axis=0)。如果 shape[i] > x.shape[i],则第 i 维度用零填充。如果 shape[i] < x.shape[i],则第 i 维度截断为长度 shape[i]。如果 shape 的任何元素为-1,则使用 x 的相应维度大小。

axes 整数或整数数组或 None,可选。

计算 DCT 的轴。默认为所有轴。

norm {None, ‘ortho’},可选。

归一化模式(见注释)。默认为 None。

overwrite_x 布尔值,可选。

如果为 True,则 x 的内容可以被破坏;默认为 False。

返回:

y 实数的 ndarray。

变换后的输入数组。

另见

idstn

多维逆正弦变换。

注释

有关 DST 类型和归一化模式的详细信息以及参考资料,请参阅 dst

示例

>>> import numpy as np
>>> from scipy.fftpack import dstn, idstn
>>> rng = np.random.default_rng()
>>> y = rng.standard_normal((16, 16))
>>> np.allclose(y, idstn(dstn(y, norm='ortho'), norm='ortho'))
True 

scipy.fftpack.idstn

原文:docs.scipy.org/doc/scipy-1.12.0/reference/generated/scipy.fftpack.idstn.html#scipy.fftpack.idstn

scipy.fftpack.idstn(x, type=2, shape=None, axes=None, norm=None, overwrite_x=False)

返回沿指定轴的多维离散正弦变换。

参数:

xarray_like

输入数组。

type,可选

DST 的类型(参见注释)。默认类型为 2。

shapeint 或整数数组或 None,可选

结果的形状。如果shapeaxes(见下文)都为 None,则shapex.shape;如果shape为 None 但axes不为 None,则shapenumpy.take(x.shape, axes, axis=0)。如果shape[i] > x.shape[i],则第 i 维用零填充。如果shape[i] < x.shape[i],则第 i 维截断为长度shape[i]。如果shape的任何元素为-1,则使用x的相应维度的大小。

axesint 或整数数组或 None,可选

计算 IDST 的轴。默认为所有轴。

norm,可选

规范化模式(参见注释)。默认为 None。

overwrite_xbool,可选

如果为 True,则x的内容可以被销毁;默认为 False。

返回:

y实数的 ndarray

转换后的输入数组。

参见

dstn

多维度 DST

注释

有关 IDST 类型和规范化模式的详细信息以及参考文献,请参见idst

示例

>>> import numpy as np
>>> from scipy.fftpack import dstn, idstn
>>> rng = np.random.default_rng()
>>> y = rng.standard_normal((16, 16))
>>> np.allclose(y, idstn(dstn(y, norm='ortho'), norm='ortho'))
True 

scipy.fftpack.diff

原文:docs.scipy.org/doc/scipy-1.12.0/reference/generated/scipy.fftpack.diff.html#scipy.fftpack.diff

scipy.fftpack.diff(x, order=1, period=None, _cache={})

返回周期序列 x 的第 k 阶导数(或积分)。

如果 x_j 和 y_j 分别是周期函数 x 和 y 的傅里叶系数,则:

y_j = pow(sqrt(-1)*j*2*pi/period, order) * x_j
y_0 = 0 if order is not 0. 

参数:

xarray_like

输入数组。

orderint,可选

差分的阶数。默认阶数为 1。如果阶数为负,则在假设 x_0 == 0 的情况下进行积分。

periodfloat,可选

序列的假设周期。默认为 2*pi

注:

如果 sum(x, axis=0) = 0,那么 diff(diff(x, k), -k) == x(在数值精度内)。

对于奇数阶和偶数 len(x),将采用 Nyquist 模式为零。

scipy.fftpack.tilbert

原文:docs.scipy.org/doc/scipy-1.12.0/reference/generated/scipy.fftpack.tilbert.html#scipy.fftpack.tilbert

scipy.fftpack.tilbert(x, h, period=None, _cache={})

返回周期序列 x 的 h-Tilbert 变换。

如果 x_j 和 y_j 是周期函数 x 和 y 的 Fourier 系数,则:

y_j = sqrt(-1)*coth(j*h*2*pi/period) * x_j
y_0 = 0 

参数:

xarray_like

要转换的输入数组。

hfloat

定义 Tilbert 变换的参数。

periodfloat, optional

序列的假定周期。默认周期为 2*pi

返回:

tilbertndarray

变换的结果。

注意

如果 sum(x, axis=0) == 0 并且 n = len(x) 是奇数,则 tilbert(itilbert(x)) == x

如果 2 * pi * h / period 大约为 10 或更大,则数值上 tilbert == hilbert(理论上 oo-Tilbert == Hilbert)。

对于偶数长度的 x,取 x 的奈奎斯特模为零。

scipy.fftpack.itilbert

原文:docs.scipy.org/doc/scipy-1.12.0/reference/generated/scipy.fftpack.itilbert.html#scipy.fftpack.itilbert

scipy.fftpack.itilbert(x, h, period=None, _cache={})

返回周期序列 x 的逆 h-Tilbert 变换。

如果 x_jy_j 是周期函数 x 和 y 的傅里叶系数,则:

y_j = -sqrt(-1)*tanh(j*h*2*pi/period) * x_j
y_0 = 0 

更多详细信息,请参见 tilbert

scipy.fftpack.hilbert

原文:docs.scipy.org/doc/scipy-1.12.0/reference/generated/scipy.fftpack.hilbert.html#scipy.fftpack.hilbert

scipy.fftpack.hilbert(x, _cache={})

返回周期序列 x 的希尔伯特变换。

如果 x_j 和 y_j 分别是周期函数 x 和 y 的傅里叶系数,则:

y_j = sqrt(-1)*sign(j) * x_j
y_0 = 0 

参数:

x 数组样式

输入数组,应该是周期性的。

_cache 字典,可选

包含用于卷积的核的字典。

返回:

y 数组

转换后的输入。

另见

scipy.signal.hilbert

使用希尔伯特变换计算分析信号。

注意事项

如果 sum(x, axis=0) == 0,那么 hilbert(ihilbert(x)) == x

对于偶数长度的 x,采用 Nyquist 模式将 x 的值设为零。

返回变换的符号没有一个常见于希尔伯特变换定义中的 -1 因子。还要注意,scipy.signal.hilbert 比这个函数多了一个 -1 因子。

scipy.fftpack.ihilbert

原文:docs.scipy.org/doc/scipy-1.12.0/reference/generated/scipy.fftpack.ihilbert.html#scipy.fftpack.ihilbert

scipy.fftpack.ihilbert(x)

返回周期序列 x 的逆 Hilbert 变换。

如果 x_jy_j 分别是周期函数 x 和 y 的 Fourier 系数,则:

y_j = -sqrt(-1)*sign(j) * x_j
y_0 = 0 

scipy.fftpack.cs_diff

原文:docs.scipy.org/doc/scipy-1.12.0/reference/generated/scipy.fftpack.cs_diff.html#scipy.fftpack.cs_diff

scipy.fftpack.cs_diff(x, a, b, period=None, _cache={})

返回周期序列的(a, b)-双曲/双曲伪导数。

如果x_jy_j是周期函数xy的傅立叶系数,则:

y_j = -sqrt(-1)*cosh(j*a*2*pi/period)/sinh(j*b*2*pi/period) * x_j
y_0 = 0 

参数:

xarray_like

要从中进行伪导数操作的数组。

a, bfloat

定义双曲/双曲伪微分算子的参数。

periodfloat, optional

序列的周期。默认周期为2*pi

返回:

cs_diffndarray

周期序列x的伪导数。

注意

对于偶数长度的xx的奈奎斯特模式被视为零。

scipy.fftpack.sc_diff

原文:docs.scipy.org/doc/scipy-1.12.0/reference/generated/scipy.fftpack.sc_diff.html#scipy.fftpack.sc_diff

scipy.fftpack.sc_diff(x, a, b, period=None, _cache={})

返回周期序列 x 的 (a,b)-双曲正弦/余弦伪导数。

如果 x_j 和 y_j 是周期函数 x 和 y 的傅里叶系数,则:

y_j = sqrt(-1)*sinh(j*a*2*pi/period)/cosh(j*b*2*pi/period) * x_j
y_0 = 0 

参数:

x类似数组

输入数组。

a,b浮点数

定义双曲正弦/余弦伪微分算子的参数。

period浮点数,可选

序列 x 的周期。默认为 2*pi。

注意

sc_diff(cs_diff(x,a,b),b,a) == x 对于偶数长度的 x,将其奈奎斯特模式视为零。

scipy.fftpack.ss_diff

原文:docs.scipy.org/doc/scipy-1.12.0/reference/generated/scipy.fftpack.ss_diff.html#scipy.fftpack.ss_diff

scipy.fftpack.ss_diff(x, a, b, period=None, _cache={})

返回(a,b)-sinh/sinh 周期序列 x 的伪导数。

如果 x 和 y 的傅里叶系数分别为 x_j 和 y_j,则:

y_j = sinh(j*a*2*pi/period)/sinh(j*b*2*pi/period) * x_j
y_0 = a/b * x_0 

参数:

x array_like

从中进行伪导数操作的数组。

a,b

定义 sinh/sinh 伪微分算子的参数。

period 浮点数,可选

序列 x 的周期。默认为 2*pi

注意事项

ss_diff(ss_diff(x,a,b),b,a) == x

scipy.fftpack.cc_diff

原文:docs.scipy.org/doc/scipy-1.12.0/reference/generated/scipy.fftpack.cc_diff.html#scipy.fftpack.cc_diff

scipy.fftpack.cc_diff(x, a, b, period=None, _cache={})

返回周期序列的 (a,b)-cosh/cosh 伪导数。

如果 x 和 y 的 Fourier 系数分别是周期函数 x 和 y 的 Fourier 系数 x_j 和 y_j,则:

y_j = cosh(j*a*2*pi/period)/cosh(j*b*2*pi/period) * x_j 

参数:

x array_like

要从伪导数中获取的数组。

a,b float

定义 sinh/sinh 伪微分算子的参数。

period float,可选

序列 x 的周期。默认为 2*pi

返回:

cc_diff ndarray

周期序列 x 的伪导数。

注意

cc_diff(cc_diff(x,a,b),b,a) == x

scipy.fftpack.shift

原文:docs.scipy.org/doc/scipy-1.12.0/reference/generated/scipy.fftpack.shift.html#scipy.fftpack.shift

scipy.fftpack.shift(x, a, period=None, _cache={})

将周期序列 x 平移 a:y(u) = x(u+a)。

如果 x_j 和 y_j 是周期函数 x 和 y 的 Fourier 系数,则:

y_j = exp(j*a*2*pi/period*sqrt(-1)) * x_f 

参数:

xarray_like

从中取伪导数的数组。

afloat

定义双曲正弦/双曲正弦伪微分的参数。

periodfloat, 可选

序列 x 和 y 的周期。默认周期是2*pi

scipy.fftpack.fftshift

原文:docs.scipy.org/doc/scipy-1.12.0/reference/generated/scipy.fftpack.fftshift.html#scipy.fftpack.fftshift

scipy.fftpack.fftshift(x, axes=None)

将零频率分量移动到频谱的中心。

此函数对所有列出的轴交换半空间(默认为所有)。注意,仅当len(x)为偶数时,y[0]才是奈奎斯特分量。

参数:

xarray_like

输入数组。

int 或者形状元组,可选

进行移动的轴。默认为 None,表示移动所有轴。

返回:

yndarray

移动后的数组。

另请参见

ifftshift

fftshift的反函数。

示例

>>> freqs = np.fft.fftfreq(10, 0.1)
>>> freqs
array([ 0.,  1.,  2., ..., -3., -2., -1.])
>>> np.fft.fftshift(freqs)
array([-5., -4., -3., -2., -1.,  0.,  1.,  2.,  3.,  4.]) 

仅沿第二轴移动零频率分量:

>>> freqs = np.fft.fftfreq(9, d=1./9).reshape(3, 3)
>>> freqs
array([[ 0.,  1.,  2.],
 [ 3.,  4., -4.],
 [-3., -2., -1.]])
>>> np.fft.fftshift(freqs, axes=(1,))
array([[ 2.,  0.,  1.],
 [-4.,  3.,  4.],
 [-1., -3., -2.]]) 

scipy.fftpack.ifftshift

原文:docs.scipy.org/doc/scipy-1.12.0/reference/generated/scipy.fftpack.ifftshift.html#scipy.fftpack.ifftshift

scipy.fftpack.ifftshift(x, axes=None)

fftshift 的反函数。尽管对于偶数长度的 x 相同,但对于奇数长度的 x,这两个函数相差一个样本。

参数:

x:array_like

输入数组。

axes:int 或者形状元组,可选

计算的轴。默认为 None,表示所有轴都移动。

返回值:

y:ndarray

移动后的数组。

参见

fftshift

将零频率分量移动到频谱中心。

示例

>>> freqs = np.fft.fftfreq(9, d=1./9).reshape(3, 3)
>>> freqs
array([[ 0.,  1.,  2.],
 [ 3.,  4., -4.],
 [-3., -2., -1.]])
>>> np.fft.ifftshift(np.fft.fftshift(freqs))
array([[ 0.,  1.,  2.],
 [ 3.,  4., -4.],
 [-3., -2., -1.]]) 

scipy.fftpack.fftshift

原文:docs.scipy.org/doc/scipy-1.12.0/reference/generated/scipy.fftpack.fftshift.html#scipy.fftpack.fftshift

scipy.fftpack.fftshift(x, axes=None)

将零频率分量移动到频谱的中心。

此函数交换所有列出的轴的半空间(默认为所有)。请注意,仅当 len(x) 为偶数时,y[0] 才是奈奎斯特分量。

参数:

xarray_like

输入数组。

axesint 或者形状元组,可选

要移动的轴。默认为 None,表示移动所有轴。

返回:

yndarray

移位数组。

另请参见

ifftshift

fftshift 的逆操作。

示例

>>> freqs = np.fft.fftfreq(10, 0.1)
>>> freqs
array([ 0.,  1.,  2., ..., -3., -2., -1.])
>>> np.fft.fftshift(freqs)
array([-5., -4., -3., -2., -1.,  0.,  1.,  2.,  3.,  4.]) 

仅沿第二轴移动零频率分量:

>>> freqs = np.fft.fftfreq(9, d=1./9).reshape(3, 3)
>>> freqs
array([[ 0.,  1.,  2.],
 [ 3.,  4., -4.],
 [-3., -2., -1.]])
>>> np.fft.fftshift(freqs, axes=(1,))
array([[ 2.,  0.,  1.],
 [-4.,  3.,  4.],
 [-1., -3., -2.]]) 

scipy.fftpack.fftfreq

原文:docs.scipy.org/doc/scipy-1.12.0/reference/generated/scipy.fftpack.fftfreq.html#scipy.fftpack.fftfreq

scipy.fftpack.fftfreq(n, d=1.0)

返回离散傅立叶变换的样本频率。

返回的浮点数数组 f 包含频率频段的中心,单位为每个样本间隔的循环数(从起始处开始)。例如,如果样本间隔以秒为单位,则频率单位为每秒循环数。

给定窗口长度 n 和样本间隔 d

f = [0, 1, ...,   n/2-1,     -n/2, ..., -1] / (d*n)   if n is even
f = [0, 1, ..., (n-1)/2, -(n-1)/2, ..., -1] / (d*n)   if n is odd 

参数:

nint

窗口长度。

dscalar, optional

样本间隔(采样率的倒数)。默认为 1。

返回:

fndarray

包含样本频率的长度为 n 的数组。

Examples

>>> signal = np.array([-2, 8, 6, 4, 1, 0, 3, 5], dtype=float)
>>> fourier = np.fft.fft(signal)
>>> n = signal.size
>>> timestep = 0.1
>>> freq = np.fft.fftfreq(n, d=timestep)
>>> freq
array([ 0\.  ,  1.25,  2.5 , ..., -3.75, -2.5 , -1.25]) 

scipy.fftpack.rfftfreq

原文:docs.scipy.org/doc/scipy-1.12.0/reference/generated/scipy.fftpack.rfftfreq.html#scipy.fftpack.rfftfreq

scipy.fftpack.rfftfreq(n, d=1.0)

DFT 样本频率(用于 rfft, irfft 的用法)。

返回的浮点数组包含频率分量(从零开始)在单位内的周期数,给定窗口长度 n 和采样间距 d

f = [0,1,1,2,2,...,n/2-1,n/2-1,n/2]/(d*n)   if n is even
f = [0,1,1,2,2,...,n/2-1,n/2-1,n/2,n/2]/(d*n)   if n is odd 

参数:

n 整数

窗口长度。

d 标量,可选

采样间距。默认为 1。

返回:

out ndarray

长度为 n 的数组,包含样本频率。

示例

>>> import numpy as np
>>> from scipy import fftpack
>>> sig = np.array([-2, 8, 6, 4, 1, 0, 3, 5], dtype=float)
>>> sig_fft = fftpack.rfft(sig)
>>> n = sig_fft.size
>>> timestep = 0.1
>>> freq = fftpack.rfftfreq(n, d=timestep)
>>> freq
array([ 0\.  ,  1.25,  1.25,  2.5 ,  2.5 ,  3.75,  3.75,  5\.  ]) 

scipy.fftpack.next_fast_len

原文:docs.scipy.org/doc/scipy-1.12.0/reference/generated/scipy.fftpack.next_fast_len.html#scipy.fftpack.next_fast_len

scipy.fftpack.next_fast_len(target)

找到输入数据的下一个快速大小,用于fft的零填充等。

SciPy 的 FFTPACK 具有基数{2, 3, 4, 5}的高效函数,因此返回大于或等于目标的下一个 2、3 和 5 的素因子的合成数。(这些也被称为 5-光滑数、正则数或 Hamming 数。)

参数:

target整数

开始搜索的长度。必须是正整数。

返回:

out整数

大于或等于目标的第一个 5-光滑数。

注意

自版本 0.18.0 起新增。

示例

在特定机器上,素数长度的 FFT 花费 133 毫秒:

>>> from scipy import fftpack
>>> import numpy as np
>>> rng = np.random.default_rng()
>>> min_len = 10007  # prime length is worst case for speed
>>> a = rng.standard_normal(min_len)
>>> b = fftpack.fft(a) 

零填充到下一个 5-光滑长度可将计算时间减少到 211 微秒,加速 630 倍:

>>> fftpack.next_fast_len(min_len)
10125
>>> b = fftpack.fft(a, 10125) 

将输入舍入到下一个 2 的幂次方并不是最优的,计算需要 367 微秒,比 5-光滑尺寸长 1.7 倍:

>>> b = fftpack.fft(a, 16384) 

scipy.fftpack.fft

原文:docs.scipy.org/doc/scipy-1.12.0/reference/generated/scipy.fftpack.fft.html#scipy.fftpack.fft

scipy.fftpack.fft(x, n=None, axis=-1, overwrite_x=False)

返回实序列或复序列的离散傅里叶变换。

返回的复数数组包含y(0), y(1),..., y(n-1),其中

y(j) = (x * exp(-2*pi*sqrt(-1)*j*np.arange(n)/n)).sum().

参数:

x类数组

要傅里叶变换的数组。

n整数,可选

傅里叶变换的长度。如果n < x.shape[axis],则截断x。如果n > x.shape[axis],则对x进行零填充。默认结果为n = x.shape[axis]

axis整数,可选

计算 FFT 的轴;默认值为最后一个轴(即,axis=-1)。

overwrite_x布尔,可选

如果为 True,则x的内容可以被破坏;默认值为 False。

返回:

z复数 ndarray

元素为:

[y(0),y(1),..,y(n/2),y(1-n/2),...,y(-1)]        if n is even
[y(0),y(1),..,y((n-1)/2),y(-(n-1)/2),...,y(-1)]  if n is odd 

其中:

y(j) = sum[k=0..n-1] x[k] * exp(-sqrt(-1)*j*k* 2*pi/n), j = 0..n-1 

另请参见

ifft

逆 FFT

rfft

实序列的 FFT

注意事项

结果的打包是“标准”的:如果A = fft(a, n),那么A[0]包含零频率项,A[1:n/2]包含正频率项,A[n/2:]按照递减负频率的顺序包含负频率项。因此,对于 8 点变换,结果的频率为[0, 1, 2, 3, -4, -3, -2, -1]。要重新排列 fft 输出以使零频率分量居中,如[-4, -3, -2, -1, 0, 1, 2, 3],请使用fftshift

实现了单精度和双精度例程。将半精度输入转换为单精度。非浮点输入将转换为双精度。不支持长双精度输入。

n是 2 的幂时,此函数最有效,当n是素数时,效率最低。

请注意,如果x是实值,则A[j] == A[n-j].conjugate()。如果x是实值且n是偶数,则A[n/2]是实数。

如果x的数据类型是实数,则自动使用“实数 FFT”算法,其大致减半计算时间。为了进一步提高效率,使用rfft,它执行相同的计算,但仅输出对称频谱的一半。如果数据既是实数又是对称的,则dct可以再次通过从信号的一半生成频谱的一半来将效率加倍。

示例

>>> import numpy as np
>>> from scipy.fftpack import fft, ifft
>>> x = np.arange(5)
>>> np.allclose(fft(ifft(x)), x, atol=1e-15)  # within numerical accuracy.
True 

scipy.fftpack.convolve.convolve

原文:docs.scipy.org/doc/scipy-1.12.0/reference/generated/scipy.fftpack.convolve.convolve.html#scipy.fftpack.convolve.convolve

scipy.fftpack.convolve.convolve(x, omega[, swap_real_imag, overwrite_x])

convolve的包装器。

参数:

x输入秩为 1 的数组(‘d’),边界为(n)

omega输入秩为 1 的数组(‘d’),边界为(n)

返回:

y秩为 1 的数组(‘d’),边界为(n),并使用 x 存储

其他参数:

overwrite_x输入整数,可选

默认:0

swap_real_imag输入整数,可选

默认:0

scipy.fftpack.convolve.convolve_z

原文:docs.scipy.org/doc/scipy-1.12.0/reference/generated/scipy.fftpack.convolve.convolve_z.html#scipy.fftpack.convolve.convolve_z

scipy.fftpack.convolve.convolve_z(x, omega_real, omega_imag[, overwrite_x])

convolve_z 的包装器。

参数:

x 输入秩-1 数组(‘d’),范围为(n)

omega_real 输入秩-1 数组(‘d’),范围为(n)

omega_imag 输入秩-1 数组(‘d’),范围为(n)

返回:

y 秩-1 数组(‘d’),范围为(n),并且具有 x 存储

其他参数:

overwrite_x 输入 int,可选

默认:0

scipy.fftpack.convolve.init_convolution_kernel

原文:docs.scipy.org/doc/scipy-1.12.0/reference/generated/scipy.fftpack.convolve.init_convolution_kernel.html#scipy.fftpack.convolve.init_convolution_kernel

scipy.fftpack.convolve.init_convolution_kernel(n, kernel_func[, d, zero_nyquist, kernel_func_extra_args])

init_convolution_kernel的包装器。

参数:

n输入整数

kernel_func回调函数

返回:

omega秩-1 数组(‘d’),边界为(n)

其他参数:

d输入整数,可选

默认值:0

kernel_func_extra_args输入元组,可选

默认值:()

zero_nyquist输入整数,可选

默认值:d%2

注释

回调函数:

def kernel_func(k): return kernel_func
Required arguments:
  k : input int
Return objects:
  kernel_func : float 

scipy.fftpack.convolve.destroy_convolve_cache

原文:docs.scipy.org/doc/scipy-1.12.0/reference/generated/scipy.fftpack.convolve.destroy_convolve_cache.html#scipy.fftpack.convolve.destroy_convolve_cache

scipy.fftpack.convolve.destroy_convolve_cache()

函数积分和常微分方程 (scipy.integrate)

原文:docs.scipy.org/doc/scipy-1.12.0/reference/integrate.html

函数积分,给定函数对象

quad(func, a, b[, args, full_output, ...]) 计算定积分。
quad_vec(f, a, b[, epsabs, epsrel, norm, ...]) 向量值函数的自适应积分。
dblquad(func, a, b, gfun, hfun[, args, ...]) 计算二重积分。
tplquad(func, a, b, gfun, hfun, qfun, rfun) 计算三重(定积分)。
nquad(func, ranges[, args, opts, full_output]) 多变量积分。
fixed_quad(func, a, b[, args, n]) 使用固定阶数的高斯积分计算定积分。
quadrature(func, a, b[, args, tol, rtol, ...]) 使用固定容差的高斯积分计算定积分。
romberg(function, a, b[, args, tol, rtol, ...]) 对可调用函数或方法进行龙贝格积分。
newton_cotes(rn[, equal]) 返回牛顿-科特斯积分的权重和误差系数。
qmc_quad(func, a, b, *[, n_estimates, ...]) 使用准蒙特卡洛积分法计算 N 维积分。
IntegrationWarning 关于积分过程中问题的警告。
AccuracyWarning

给定固定样本的函数积分

trapezoid(y[, x, dx, axis]) 使用复合梯形法则沿给定轴积分。
cumulative_trapezoid(y[, x, dx, axis, initial]) 使用复合梯形法累积积分 y(x)。
simpson(y, *[, x, dx, axis, even]) 使用给定轴上的样本和复合 Simpson 法积分 y(x)。
cumulative_simpson(y, *[, x, dx, axis, initial]) 使用复合 Simpson's 1/3 法累积积分 y(x)。
romb(y[, dx, axis, show]) 使用函数样本的 Romberg 积分。

另见

scipy.special 用于正交多项式(特殊函数)的高斯积分根和其他权重因子和区域。

解决 ODE 系统的初值问题

这些求解器被实现为各自的类,可以直接使用(低级用法)或通过便捷函数使用。

solve_ivp(fun, t_span, y0[, method, t_eval, ...]) 解决 ODE 系统的初值问题。
RK23(fun, t0, y0, t_bound[, max_step, rtol, ...]) 3(2)阶显式 Runge-Kutta 方法。
RK45(fun, t0, y0, t_bound[, max_step, rtol, ...]) 5(4)阶显式 Runge-Kutta 方法。
DOP853(fun, t0, y0, t_bound[, max_step, ...]) 8 阶显式 Runge-Kutta 方法。
Radau(fun, t0, y0, t_bound[, max_step, ...]) Radau IIA 家族的隐式 Runge-Kutta 方法,5 阶。
BDF(fun, t0, y0, t_bound[, max_step, rtol, ...]) 基于后向差分公式的隐式方法。
LSODA(fun, t0, y0, t_bound[, first_step, ...]) 具有自动刚度检测和切换的 Adams/BDF 方法。
OdeSolver(fun, t0, y0, t_bound, vectorized) ODE 求解器的基类。
DenseOutput 用于 ODE 求解器在步长上的局部插值的基类。
OdeSolution 连续的 ODE 解。

旧 API

这些是早期为 SciPy 开发的例程。它们封装了用 Fortran 实现的旧求解器(主要是 ODEPACK)。虽然它们的接口并不特别方便,与新 API 相比某些功能也不完整,但这些求解器本身质量良好且作为编译后的 Fortran 代码运行速度快。在某些情况下,使用这个旧 API 可能是值得的。

odeint 积分一组常微分方程。
ode 一个通用的数值积分器接口类。
complex_ode 复杂系统的 ODE 包装器。
ODEintWarning 在执行 odeint 过程中引发的警告。

解决常微分方程组的边界值问题。

solve_bvp 解决常微分方程组的边界值问题。

scipy.integrate.quad

原文:docs.scipy.org/doc/scipy-1.12.0/reference/generated/scipy.integrate.quad.html#scipy.integrate.quad

scipy.integrate.quad(func, a, b, args=(), full_output=0, epsabs=1.49e-08, epsrel=1.49e-08, limit=50, points=None, weight=None, wvar=None, wopts=None, maxp1=50, limlst=50, complex_func=False)

计算定积分。

使用 Fortran 库 QUADPACK 中的技术从ab(可能是无限区间)积分func

参数:

func

用于积分的 Python 函数或方法。如果func接受多个参数,则沿着与第一个参数对应的轴积分。

如果用户希望改进积分性能,则f可以是具有以下签名之一的scipy.LowLevelCallable

double func(double x)
double func(double x, void *user_data)
double func(int n, double *xx)
double func(int n, double *xx, void *user_data) 

user_data是包含在scipy.LowLevelCallable中的数据。在带有xx的调用形式中,nxx数组的长度,其中包含xx[0] == x,其余项目是 quad 函数的args参数中包含的数字。

此外,某些 ctypes 调用签名支持向后兼容性,但不应在新代码中使用。

a浮点数

积分的下限(使用-numpy.inf 表示-无穷大)。

b浮点数

积分的上限(使用 numpy.inf 表示+无穷大)。

args元组,可选

额外传递给func的参数。

full_output整数,可选

非零以返回积分信息的字典。如果非零,则还抑制警告消息并将消息追加到输出元组中。

complex_func布尔值,可选

指示函数(func)返回类型是否为实数(complex_func=False:默认)或复数(complex_func=True)。在两种情况下,函数的参数是实数。如果full_output也非零,则实部和虚部的infodictmessageexplain以“real output”和“imag output”为键返回到字典中。

返回:

y浮点数

ab的函数func的积分。

abserr浮点数

结果的绝对误差估计。

infodict字典

包含附加信息的字典。

消息

收敛消息。

解释

仅在具有“cos”或“sin”加权和无限积分限制时附加,它包含 infodict['ierlst']中代码的解释。

其他参数:

epsabs浮点数或整数,可选

绝对误差容限。默认为 1.49e-8。quad试图获得abs(i-result) <= max(epsabs, epsrel*abs(i))的精度,其中i = funcab的积分,而result是数值近似值。见下文的epsrel

epsrel浮点数或整数,可选

相对误差容限。默认为 1.49e-8。如果epsabs <= 0epsrel必须大于 5e-29 和50 * (machine epsilon)。见上述epsabs

limit浮点数或整数,可选

自适应算法中使用的子区间数量的上限。

points(sequence of floats,ints), optional

有界积分区间中可能发生积分被积函数的局部困难(例如奇点、不连续点)的断点序列。序列不必排序。请注意,此选项不能与 weight 结合使用。

weightfloat or int, optional

指示加权函数的字符串。有关此及其余参数的详细说明,请参阅下文。

wvaroptional

变量,用于加权函数。

woptsoptional

重复使用切比雪夫矩的可选输入。

maxp1float or int, optional

切比雪夫矩的数量上限。

limlstint, optional

循环数量的上限(>=3)适用于正弦加权和无限端点。

另请参阅

dblquad

双重积分

tplquad

三重积分

nquad

n 维积分(递归使用 quad

fixed_quad

固定阶数的高斯积分

quadrature

自适应高斯积分

odeint

ODE 积分器

ode

ODE 积分器

simpson

采样数据的积分器

romb

采样数据的积分器

scipy.special

用于正交多项式的系数和根

注意事项

积分必须收敛以获得有效结果;不保证发散积分的行为。

quad() 输入和输出的额外信息

如果 full_output 非零,则第三个输出参数(infodict)是一个具有如下表格条目的字典。对于无限限制,范围转换为 (0,1),并给出了相对于此转换范围的可选输出。令 M 为输入参数限制,K 为 infodict[‘last’]。条目如下:

‘neval’

函数评估的数量。

‘last’

分割过程中产生的子区间数量 K。

‘alist’

长度为 M 的秩-1 数组,其前 K 个元素是积分范围内分区的左端点。

‘blist’

长度为 M 的秩-1 数组,其前 K 个元素是子区间的右端点。

‘rlist’

一个长度为 M 的一维数组,其前 K 个元素是子区间上的积分近似值。

‘elist’

一个长度为 M 的一维数组,其前 K 个元素是子区间上的绝对误差估计的模。

‘iord’

一个长度为 M 的一维整数数组,其前 L 个元素是子区间上的误差估计的指针,如果K<=M/2+2,则L=K,否则L=M+1-K。设 I 为序列infodict['iord'],E 为序列infodict['elist'],则E[I[1]], ..., E[I[L]]形成一个递减序列。

如果提供了输入参数 points(即它不是 None),则将以下额外输出放置在输出字典中。假设 points 序列的长度为 P。

‘pts’

一个长度为 P+2 的一维数组,按升序给出积分限和区间的断点。这是一个数组,提供将发生积分的子区间。

‘level’

一个长度为 M 的一维整数数组(即 limit),包含子区间的分割级别,即如果(aa,bb)是pts[1], pts[2]之间的子区间,其中pts[0]pts[2]infodict['pts']的相邻元素,则(aa,bb)的级别为 l,如果|bb-aa| = |pts[2]-pts[1]| * 2**(-l)

‘ndin’

一个长度为 P+2 的一维整数数组。在第一次积分后,一些区间的误差估计可能会被人为增加,以推动它们的分割。这个数组在对应于发生这种情况的子区间的槽中有 1。

加权积分

输入变量weightwvar用于通过一组选择的函数对积分被加权。使用这些加权函数计算积分的不同方法,并且这些不支持指定断点。weight 的可能值及其对应的加权函数如下。

weight 使用的加权函数 wvar
‘cos’ cos(w*x) wvar = w
‘sin’ sin(w*x) wvar = w
‘alg’ g(x) = ((x-a)alpha)*((b-x)beta) wvar = (alpha, beta)
‘alg-loga’ g(x)*log(x-a) wvar = (alpha, beta)
‘alg-logb’ g(x)*log(b-x) wvar = (alpha, beta)
‘alg-log’ g(x)log(x-a)log(b-x) wvar = (alpha, beta)
‘cauchy’ 1/(x-c) wvar = c

这些表达式中,a 和 b 是积分限。

对于‘cos’和‘sin’加权,提供了额外的输入和输出。

对于有限的积分限,使用 Clenshaw-Curtis 方法执行积分,该方法使用切比雪夫矩。对于重复计算,这些矩保存在输出字典中:

‘momcom’

已计算的切比雪夫矩数的最大级别,即如果M_cinfodict['momcom'],则已对长度为|b-a| * 2**(-l)的区间(其中l=0,1,...,M_c)进行了计算。

| ‘nnlog’

一个长度为 M(=limit)的秩为 1 的整数数组,包含子区间的分割级别,即,如果这个数组的一个元素等于 l,那么相应的子区间就是|b-a|* 2**(-l)

‘chebmo’

一个形状为(25, maxp1)的秩为 2 的数组,包含计算得到的切比雪夫矩。可以通过将此数组作为序列 wopts 的第二个元素并将 infodict['momcom']作为第一个元素,将这些传递到相同区间的积分。

如果一个积分限制为无穷大,则计算傅里叶积分(假设 w neq 0)。如果 full_output 为 1 且遇到数值错误,则除了附加到输出元组的错误消息之外,还会附加一个字典到输出元组,该字典将数组info['ierlst']中的错误代码翻译为英文消息。输出信息字典包含以下条目,而不是‘last’,‘alist’,‘blist’,‘rlist’和‘elist’:

‘lst’

积分所需的子区间数目(称之为K_f)。

‘rslst’

一个长度为 M_f=limlst 的秩为 1 的数组,其前K_f个元素包含区间(a+(k-1)c, a+kc)上的积分贡献,其中c = (2*floor(|w|) + 1) * pi / |w|k=1,2,...,K_f

‘erlst’

一个长度为M_f的秩为 1 的数组,包含与infodict['rslist']中相同位置的区间对应的误差估计。

‘ierlst’

一个长度为M_f的秩为 1 的整数数组,包含与infodict['rslist']中相同位置的区间对应的错误标志。查看输出元组中的解释字典(最后一个条目)以获取代码含义。

QUADPACK 级别例程的详细信息

quad调用来自 FORTRAN 库 QUADPACK 的例程。本节提供了每个例程被调用的条件以及每个例程的简短描述。调用的例程取决于weightpoints和积分限制ab

QUADPACK 例程 weight points 无限边界
qagse
qagie
qagpe
qawoe ‘sin’, ‘cos’
qawfe ‘sin’, ‘cos’ 要么a要么b
qawse ‘alg*’
qawce ‘cauchy’

以下从[1]提供了每个例程的简短描述。

qagse

是一种基于全局自适应区间细分和外推的积分器,将消除几种类型的被积函数奇点的影响。

qagie

处理无限区间上的积分。无限范围映射到有限区间,随后采用与QAGS相同的策略。

qagpe

具有与 QAGS 相同目的的服务,但还允许用户提供关于麻烦点位置和类型的明确信息,即内部奇点,间断点和被积函数的其他困难。

qawoe

是对在有限区间([a,b])上评估(\int^b_a \cos(\omega x)f(x)dx)或(\int^b_a \sin(\omega x)f(x)dx)的积分器,其中用户指定(\omega)和(f)。规则评估组件基于修改的 Clenshaw-Curtis 技术

自适应细分方案与外推程序结合使用,这是QAGS中的修改,允许算法处理(f(x))中的奇点。

qawfe

计算傅里叶变换(\int^\infty_a \cos(\omega x)f(x)dx)或(\int^\infty_a \sin(\omega x)f(x)dx),用户提供(\omega)和(f)。QAWO的过程应用于连续的有限区间,通过(\varepsilon)-算法对积分逼近序列进行收敛加速。

qawse

近似计算(\int^b_a w(x)f(x)dx),其中(a < b),其中(w(x) = (x-a){\alpha}(b-x)v(x)),(\alpha,\beta > -1),其中(v(x))可能是以下函数之一:(1)、(\log(x-a))、(\log(b-x))、(\log(x-a)\log(b-x))。

用户指定(\alpha)、(\beta)和函数(v)的类型。采用全局自适应细分策略,在包含ab的子区间上进行修改的 Clenshaw-Curtis 积分。

qawce

计算(\int^b_a f(x) / (x-c)dx),其中积分必须解释为柯西主值积分,用户指定(c)和(f)。采用全局自适应策略。在包含点(x = c)的那些区间上使用修改的 Clenshaw-Curtis 积分。

实变量的复函数积分

一个实变量的复值函数(f)可以写成(f = g + ih)。类似地,(f)的积分可以写成

[\int_a^b f(x) dx = \int_a^b g(x) dx + i\int_a^b h(x) dx]

假设(g)和(h)在区间([a,b])上的积分存在[2]。因此,quad通过分别积分实部和虚部来积分复值函数。

参考文献

[1]

Piessens, Robert; de Doncker-Kapenga, Elise; Überhuber, Christoph W.; Kahaner, David (1983). QUADPACK:用于自动积分的子程序包。Springer-Verlag. ISBN 978-3-540-12553-2.

[2]

McCullough, Thomas; Phillips, Keith (1973). Foundations of Analysis in the Complex Plane. Holt Rinehart Winston. ISBN 0-03-086370-8

例子

计算(\int⁴_0 x² dx)并与解析结果比较

>>> from scipy import integrate
>>> import numpy as np
>>> x2 = lambda x: x**2
>>> integrate.quad(x2, 0, 4)
(21.333333333333332, 2.3684757858670003e-13)
>>> print(4**3 / 3.)  # analytical result
21.3333333333 

计算(\int^\infty_0 e^{-x} dx)

>>> invexp = lambda x: np.exp(-x)
>>> integrate.quad(invexp, 0, np.inf)
(1.0, 5.842605999138044e-11) 

计算(\int¹_0 a x ,dx),其中(a = 1, 3)

>>> f = lambda x, a: a*x
>>> y, err = integrate.quad(f, 0, 1, args=(1,))
>>> y
0.5
>>> y, err = integrate.quad(f, 0, 1, args=(3,))
>>> y
1.5 

用 ctypes 计算(\int¹_0 x² + y² dx),其中 y 参数为 1:

testlib.c =>
    double func(int n, double args[n]){
        return args[0]*args[0] + args[1]*args[1];}
compile to library testlib.* 
from scipy import integrate
import ctypes
lib = ctypes.CDLL('/home/.../testlib.*') #use absolute path
lib.func.restype = ctypes.c_double
lib.func.argtypes = (ctypes.c_int,ctypes.c_double)
integrate.quad(lib.func,0,1,(1))
#(1.3333333333333333, 1.4802973661668752e-14)
print((1.0**3/3.0 + 1.0) - (0.0**3/3.0 + 0.0)) #Analytic result
# 1.3333333333333333 

请注意,与积分区间的尺寸相比,脉冲形状和其他尖锐特征可能无法使用这种方法正确积分。一个简化的例子是在积分边界内具有许多零值的 y 轴反射阶跃函数的积分。

>>> y = lambda x: 1 if x<=0 else 0
>>> integrate.quad(y, -1, 1)
(1.0, 1.1102230246251565e-14)
>>> integrate.quad(y, -1, 100)
(1.0000000002199108, 1.0189464580163188e-08)
>>> integrate.quad(y, -1, 10000)
(0.0, 0.0) 

scipy.integrate.quad_vec

原文:docs.scipy.org/doc/scipy-1.12.0/reference/generated/scipy.integrate.quad_vec.html#scipy.integrate.quad_vec

scipy.integrate.quad_vec(f, a, b, epsabs=1e-200, epsrel=1e-08, norm='2', cache_size=100000000.0, limit=10000, workers=1, points=None, quadrature=None, full_output=False, *, args=())

向量值函数的自适应积分。

参数:

f可调用对象

要积分的向量值函数 f(x)。

a浮点数

起点。

b浮点数

终点。

epsabs浮点数,可选

绝对容差。

epsrel浮点数,可选

相对容差。

norm,可选

用于误差估计的向量范数。

cache_size整数,可选

用于记忆化的字节数。

limit浮点数或整数,可选

自适应算法中使用的子区间数量的上限。

workers整数或类似映射的可调用对象,可选

如果workers是整数,则部分计算以并行方式划分为这么多任务(使用multiprocessing.pool.Pool)。提供-1以使用进程可用的所有核心。或者,提供一个类似映射的可调用对象,如multiprocessing.pool.Pool.map,用于并行评估人口。此评估作为workers(func, iterable)执行。

points列表,可选

附加断点列表。

quadrature,可选

在子区间上使用的积分规则。选项:‘gk21’(Gauss-Kronrod 21 点规则),‘gk15’(Gauss-Kronrod 15 点规则),‘trapezoid’(复合梯形规则)。默认值:对有限区间使用‘gk21’,对(半)无限区间使用‘gk15’。

full_output布尔型,可选

返回额外的info字典。

args元组,可选

如有需要,传递给函数的额外参数。

自 1.8.0 版本新功能。

返回值:

res

结果的估计值

err浮点数

在给定范数下结果的误差估计。

info字典

仅在full_output=True时返回。信息字典。是一个具有以下属性的对象:

成功标志布尔型

是否达到了目标精度。

status 整数

收敛的指示器,成功(0),失败(1),以及由于舍入误差而失败(2)。

neval 整数

函数评估的数量。

intervals 数组,形状(num_intervals,2)

子区间的起始点和结束点。

integrals 数组,形状(num_intervals,…)

每个区间的积分。请注意,最多记录cache_size个值,并且数组可能包含nan表示缺失项。

errors 数组,形状(num_intervals,)

每个区间的估计积分误差。

注意事项

该算法主要遵循 QUADPACK 的 DQAG*算法的实现,实现全局误差控制和自适应细分。

此处的算法与 QUADPACK 方法略有不同:

算法不是一次性将一个区间细分,而是一次性将具有最大误差的 N 个区间细分。这使得积分的(部分)并行化成为可能。

然后,不实现“下一个最大”区间优先细分的逻辑,我们依赖上述扩展来避免仅集中在“小”区间上。

Wynn epsilon 表外推法未被使用(QUADPACK 用于无限区间)。这是因为这里的算法应该适用于矢量值函数,在用户指定的范数下,而将 epsilon 算法扩展到这种情况似乎并没有得到广泛认可。对于最大范数,使用逐元素 Wynn epsilon 可能是可能的,但我们在这里没有这样做,希望 epsilon 外推法主要在特殊情况下有用。

参考文献

[1] R. Piessens, E. de Doncker, QUADPACK (1983).

例子

我们可以计算矢量值函数的积分:

>>> from scipy.integrate import quad_vec
>>> import numpy as np
>>> import matplotlib.pyplot as plt
>>> alpha = np.linspace(0.0, 2.0, num=30)
>>> f = lambda x: x**alpha
>>> x0, x1 = 0, 2
>>> y, err = quad_vec(f, x0, x1)
>>> plt.plot(alpha, y)
>>> plt.xlabel(r"$\alpha$")
>>> plt.ylabel(r"$\int_{0}^{2} x^\alpha dx$")
>>> plt.show() 

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

scipy.integrate.dblquad

原文:docs.scipy.org/doc/scipy-1.12.0/reference/generated/scipy.integrate.dblquad.html#scipy.integrate.dblquad

scipy.integrate.dblquad(func, a, b, gfun, hfun, args=(), epsabs=1.49e-08, epsrel=1.49e-08)

计算双重积分。

返回func(y, x)x = a..by = gfun(x)..hfun(x)的双(确定)积分。

参数:

func可调用。

一个 Python 函数或至少两个变量的方法:y 必须是第一个参数,x 是第二个参数。

a, b浮点数。

x 的积分限制:a < b

gfun可调用或浮点数。

y 的下边界曲线,它是一个接受单个浮点参数(x)并返回浮点结果或指示常数边界曲线的浮点数。

hfun可调用或浮点数。

y 的上边界曲线(与gfun具有相同要求)。

args序列,可选。

传递给func的额外参数。

epsabs浮点数,可选。

直接传递给内部 1-D 积分的绝对容差。默认为 1.49e-8。dblquad试图获得abs(i-result) <= max(epsabs, epsrel*abs(i))的精度,其中ifunc(y, x)gfun(x)hfun(x)的内积分,result是数值近似值。见下面的epsrel

epsrel浮点数,可选。

内部 1-D 积分的相对容差。默认为 1.49e-8。如果epsabs <= 0epsrel必须大于 5e-29 和50 * (machine epsilon)。见上面的epsabs

返回:

y浮点数。

结果积分。

abserr浮点数。

误差的估计。

另请参阅

quad

单重积分。

tplquad

三重积分。

nquad

N 维积分。

fixed_quad

固定阶高斯积分。

quadrature

自适应高斯积分。

odeint

ODE(常微分方程)积分器。

ode

ODE(常微分方程)积分器。

simpson

用于采样数据的积分器。

romb

用于采样数据的积分器。

scipy.special

用于正交多项式的系数和根。

注意:

为了获得有效的结果,积分必须收敛;对于发散的积分,行为不能保证。

QUADPACK 级别例程的详细信息

quad 调用来自 FORTRAN 库 QUADPACK 的例程。本节详细介绍了调用每个例程的条件以及每个例程的简短描述。对于每个积分级别,如果限制是有限的,则使用 qagse,如果任一限制(或两者!)是无限的,则使用 qagie。以下提供了来自 [1] 的每个例程的简短描述。

qagse

是基于全局自适应区间细分与外推结合的积分器,将消除多种类型的被积函数奇点的影响。

qagie

处理无限区间上的积分。无限范围被映射到有限区间,随后采用与 QAGS 相同的策略。

参考文献

[1]

Piessens, Robert; de Doncker-Kapenga, Elise; Überhuber, Christoph W.; Kahaner, David(1983)。QUADPACK:用于自动积分的子程序包。Springer-Verlag。ISBN 978-3-540-12553-2。

示例

计算 x * y**2 在区间 x 从 0 到 2,y 从 0 到 1 的双重积分。即 (\int^{x=2}{x=0} \int^{y=1} x y² ,dy ,dx)。

>>> import numpy as np
>>> from scipy import integrate
>>> f = lambda y, x: x*y**2
>>> integrate.dblquad(f, 0, 2, 0, 1)
 (0.6666666666666667, 7.401486830834377e-15) 

计算 (\int^{x=\pi/4}{x=0} \int^{y=\cos(x)} 1 ,dy ,dx)。

>>> f = lambda y, x: 1
>>> integrate.dblquad(f, 0, np.pi/4, np.sin, np.cos)
 (0.41421356237309503, 1.1083280054755938e-14) 

计算 (\int^{x=1}{x=0} \int^{y=2-x} a x y ,dy ,dx),其中 (a=1, 3)。

>>> f = lambda y, x, a: a*x*y
>>> integrate.dblquad(f, 0, 1, lambda x: x, lambda x: 2-x, args=(1,))
 (0.33333333333333337, 5.551115123125783e-15)
>>> integrate.dblquad(f, 0, 1, lambda x: x, lambda x: 2-x, args=(3,))
 (0.9999999999999999, 1.6653345369377348e-14) 

计算二维高斯积分,即高斯函数 (f(x,y) = e{-(x + y^{2})}) 在 ((-\infty,+\infty)) 上的积分。即计算积分 (\iint^{+\infty}_{-\infty} e{-(x + y^{2})} ,dy,dx)。

>>> f = lambda x, y: np.exp(-(x ** 2 + y ** 2))
>>> integrate.dblquad(f, -np.inf, np.inf, -np.inf, np.inf)
 (3.141592653589777, 2.5173086737433208e-08) 

scipy.integrate.tplquad

原文:docs.scipy.org/doc/scipy-1.12.0/reference/generated/scipy.integrate.tplquad.html#scipy.integrate.tplquad

scipy.integrate.tplquad(func, a, b, gfun, hfun, qfun, rfun, args=(), epsabs=1.49e-08, epsrel=1.49e-08)

计算三重(确定)积分。

返回func(z, y, x)x = a..by = gfun(x)..hfun(x),和z = qfun(x,y)..rfun(x,y)的三重积分。

参数:

func函数

一个 Python 函数或至少三个变量的方法,顺序为(z,y,x)。

a, b浮点数

x 的积分限制:a < b

gfun函数或浮点数

y 中的下边界曲线,它是一个函数,接受单个浮点参数(x)并返回浮点结果或表示常数边界曲线的浮点数。

hfun函数或浮点数

y 中的上边界曲线(与gfun要求相同)。

qfun函数或浮点数

z 中的下边界面。它必须是一个函数,接受顺序为(x,y)的两个浮点数,并返回一个浮点数或表示常数边界面的浮点数。

rfun函数或浮点数

z 中的上边界面。(与qfun要求相同。)

args元组,可选

传递给func的额外参数。

epsabs浮点数,可选

直接传递给最内层的一维积分的绝对容差。默认值为 1.49e-8。

epsrel浮点数,可选

最内层一维积分的相对容差。默认值为 1.49e-8。

返回:

y浮点数

结果积分。

abserr浮点数

误差的估计。

另请参见

quad

使用 QUADPACK 的自适应积分

quadrature

自适应高斯积分

fixed_quad

固定阶高斯积分

dblquad

双重积分

nquad

N 维积分

romb

采样数据的积分器

simpson

采样数据的积分器

ode

ODE 积分器

odeint

ODE 积分器

scipy.special

用于正交多项式的系数和根

注意事项

为了获得有效的结果,积分必须收敛;不保证发散积分的行为。

QUADPACK 级别例程的详细信息

quad 调用来自 FORTRAN 库 QUADPACK 的例程。本节提供每个例程调用条件的详细说明以及每个例程的简短描述。对于每个积分级别,如果限制是有限的,使用 qagse;如果任一限制(或两个限制!)是无限的,则使用 qagie。以下提供了来自 [1] 的每个例程的简短描述。

qagse

是一种基于全局自适应区间细分的积分器,结合外推法,可以消除多种类型的被积函数奇异性的影响。

qagie

处理对无限区间的积分。无限范围映射到有限区间,随后采用与 QAGS 相同的策略。

参考文献

[1]

Piessens, Robert; de Doncker-Kapenga, Elise; Überhuber, Christoph W.; Kahaner, David (1983). QUADPACK: A subroutine package for automatic integration. Springer-Verlag. ISBN 978-3-540-12553-2.

例子

计算三重积分 x * y * z,其中 x 范围从 1 到 2,y 范围从 2 到 3,z 范围从 0 到 1。即,(\int^{x=2}{x=1} \int^{y=3} \int^{z=1}_{z=0} x y z ,dz ,dy ,dx)。

>>> import numpy as np
>>> from scipy import integrate
>>> f = lambda z, y, x: x*y*z
>>> integrate.tplquad(f, 1, 2, 2, 3, 0, 1)
(1.8749999999999998, 3.3246447942574074e-14) 

计算 (\int^{x=1}{x=0} \int^{y=1-2x} \int^{z=1-x-2y}_{z=0} x y z ,dz ,dy ,dx)。注意:qfun/rfun 按顺序 (x, y) 接受参数,即使 f 按顺序 (z, y, x) 接受参数。

>>> f = lambda z, y, x: x*y*z
>>> integrate.tplquad(f, 0, 1, 0, lambda x: 1-2*x, 0, lambda x, y: 1-x-2*y)
(0.05416666666666668, 2.1774196738157757e-14) 

计算 (\int^{x=1}{x=0} \int^{y=1} \int^{z=1}_{z=0} a x y z ,dz ,dy ,dx) 对于 (a=1, 3)。

>>> f = lambda z, y, x, a: a*x*y*z
>>> integrate.tplquad(f, 0, 1, 0, 1, 0, 1, args=(1,))
 (0.125, 5.527033708952211e-15)
>>> integrate.tplquad(f, 0, 1, 0, 1, 0, 1, args=(3,))
 (0.375, 1.6581101126856635e-14) 

计算三维高斯积分,即高斯函数 (f(x,y,z) = e{-(x + y^{2} + z^{2})}) 在 ((-\infty,+\infty)) 上的积分。即,计算积分 (\iiint^{+\infty}_{-\infty} e{-(x + y^{2} + z^{2})} ,dz ,dy,dx)。

>>> f = lambda x, y, z: np.exp(-(x ** 2 + y ** 2 + z ** 2))
>>> integrate.tplquad(f, -np.inf, np.inf, -np.inf, np.inf, -np.inf, np.inf)
 (5.568327996830833, 4.4619078828029765e-08) 

scipy.integrate.nquad

原文:docs.scipy.org/doc/scipy-1.12.0/reference/generated/scipy.integrate.nquad.html#scipy.integrate.nquad

scipy.integrate.nquad(func, ranges, args=None, opts=None, full_output=False)

对多个变量进行积分。

包装quad以便对多个变量进行积分。各种选项允许改进不连续函数的积分,以及使用加权积分,通常更好地控制积分过程。

参数:

func {可调用对象, scipy.LowLevelCallable}

要进行积分的函数。具有x0, ... xnt0, ... tm的参数,其中积分是在x0, ... xn上进行的,这些必须是浮点数。其中t0, ... tm是通过 args 传递的额外参数。函数签名应为func(x0, x1, ..., xn, t0, t1, ..., tm)。积分是按顺序进行的。即,对x0的积分是最内层积分,而xn是最外层。

如果用户希望改进积分性能,则f可以是带有以下签名之一的scipy.LowLevelCallable

double func(int n, double *xx)
double func(int n, double *xx, void *user_data) 

其中n是变量和参数的数量。xx数组包含坐标和额外参数。user_data是包含在scipy.LowLevelCallable中的数据。

ranges 可迭代对象

ranges 的每个元素可以是 2 个数字的序列,或者是返回这样一个序列的可调用对象。ranges[0]对应于对 x0 的积分,依此类推。如果 ranges 的一个元素是可调用的,则它将使用所有可用的积分参数以及任何参数化参数进行调用。例如,如果func = f(x0, x1, x2, t0, t1),那么ranges[0]可以定义为(a, b)或者(a, b) = range0(x1, x2, t0, t1)

args 可迭代对象,可选

funcrangesopts要求的额外参数t0, ... tn

opts 可迭代对象或字典,可选

要传递给quad的选项。可以为空、字典或返回字典或函数序列。如果为空,则使用 scipy.integrate.quad 的默认选项。如果是字典,则所有积分级别使用相同的选项。如果是序列,则序列的每个元素对应于特定积分。例如,opts[0]对应于对x0的积分,依此类推。如果是可调用的,则签名必须与ranges相同。可用选项及其默认值如下:

  • epsabs = 1.49e-08
  • epsrel = 1.49e-08
  • limit = 50
  • points = None
  • weight = None
  • wvar = None
  • wopts = None

关于这些选项的更多信息,请参见quad

full_output 布尔值,可选

来自 scipy.integrate.quad 的 full_output 的部分实现。通过在调用 nquad 时设置 full_output=True 可以获取积分函数 neval 的数量。

返回:

resultfloat

积分结果。

abserrfloat

在各种积分结果的绝对误差估计的最大值。

out_dictdict,可选

包含有关积分附加信息的字典。

另请参阅

quad

1-D 数值积分

dblquad, tplquad

双重和三重积分

fixed_quad

固定阶数的高斯积分

quadrature

自适应高斯积分

注意

为了获得有效结果,积分必须收敛;对于发散的积分,结果不能保证。

QUADPACK 等级例程的详细信息

nquad 调用来自 FORTRAN 库 QUADPACK 的例程。本节提供了每个例程被调用的条件和每个例程的简短描述。所调用的例程取决于 weightpoints 和积分限 ab

QUADPACK 程序 weight points 无限界限
qagse
qagie
qagpe
qawoe ‘sin’, ‘cos’
qawfe ‘sin’, ‘cos’ ab 中的任一者
qawse ‘alg*’
qawce ‘cauchy’

以下提供了每个例程的简短描述,来源于[1]

qagse

是基于全局自适应区间分割与外推结合的积分器,它将消除几种类型积分函数奇点的影响。

qagie

处理无限区间上的积分。将无限范围映射到有限区间,随后应用与 QAGS 中相同的策略。

qagpe

与 QAGS 有相同的功能,但也允许用户提供关于麻烦点(如积分函数内部奇异性、不连续性和其他难点的抛物线的位置和类型的明确信息。

qawoe

是一个用于计算 (\int^b_a \cos(\omega x)f(x)dx) 或 (\int^b_a \sin(\omega x)f(x)dx) 的积分器,其中用户指定了 (\omega) 和 (f)。规则评估组件基于修改的 Clenshaw-Curtis 技术。

使用与 QAGS 中的修改相同的外推程序的自适应分段方案,这将消除几种类型的积分函数奇点的影响。

qawfe

计算用户提供的 (\omega) 和 (f) 的傅里叶变换 (\int^\infty_a \cos(\omega x)f(x)dx) 或 (\int^\infty_a \sin(\omega x)f(x)dx)。QAWO 过程应用于连续的有限区间,通过 (\varepsilon)-算法加速收敛到积分逼近的级数。

qawse

近似计算 (\int^b_a w(x)f(x)dx),其中 (a < b),(w(x) = (x-a){\alpha}(b-x)v(x)),(\alpha,\beta > -1),(v(x)) 可能是以下函数之一:(1),(\log(x-a)),(\log(b-x)),(\log(x-a)\log(b-x))。

用户指定 (\alpha)、(\beta) 和函数 (v) 的类型。应用全局自适应细分策略,在包含 ab 的子区间上使用改进的 Clenshaw-Curtis 积分。

qawce

计算 (\int^b_a f(x) / (x-c)dx),积分必须解释为柯西主值积分,对于用户指定的 (c) 和 (f)。采用全局自适应策略。在包含点 (x = c) 的区间上使用改进的 Clenshaw-Curtis 积分。

参考文献

[1]

Piessens, Robert; de Doncker-Kapenga, Elise; Überhuber, Christoph W.; Kahaner, David (1983). QUADPACK: 一个用于自动积分的子程序包。Springer-Verlag。ISBN 978-3-540-12553-2。

示例

计算

[\int^{1}{-0.15} \int^{0.8} \int^{1}{-1} \int^{1} f(x_0, x_1, x_2, x_3) ,dx_0 ,dx_1 ,dx_2 ,dx_3 ,]

其中

[\begin{split}f(x_0, x_1, x_2, x_3) = \begin{cases} x_0²+x_1 x_2-x_3³+ \sin{x_0}+1 & (x_0-0.2 x_3-0.5-0.25 x_1 > 0) \ x_0²+x_1 x_2-x_3³+ \sin{x_0}+0 & (x_0-0.2 x_3-0.5-0.25 x_1 \leq 0) \end{cases} .\end{split}]

>>> import numpy as np
>>> from scipy import integrate
>>> func = lambda x0,x1,x2,x3 : x0**2 + x1*x2 - x3**3 + np.sin(x0) + (
...                                 1 if (x0-.2*x3-.5-.25*x1>0) else 0)
>>> def opts0(*args, **kwargs):
...     return {'points':[0.2*args[2] + 0.5 + 0.25*args[0]]}
>>> integrate.nquad(func, [[0,1], [-1,1], [.13,.8], [-.15,1]],
...                 opts=[opts0,{},{},{}], full_output=True)
(1.5267454070738633, 2.9437360001402324e-14, {'neval': 388962}) 

计算

[\int^{t_0+t_1+1}{t_0+t_1-1} \int^{x_2+t_0² t_1³+1} \int^{t_0 x_1+t_1 x_2+1}_{t_0 x_1+t_1 x_2-1} f(x_0,x_1, x_2,t_0,t_1) ,dx_0 ,dx_1 ,dx_2,]

其中

[\begin{split}f(x_0, x_1, x_2, t_0, t_1) = \begin{cases} x_0 x_2² + \sin{x_1}+2 & (x_0+t_1 x_1-t_0 > 0) \ x_0 x_2² +\sin{x_1}+1 & (x_0+t_1 x_1-t_0 \leq 0) \end{cases}\end{split}]

和 ((t_0, t_1) = (0, 1))。

>>> def func2(x0, x1, x2, t0, t1):
...     return x0*x2**2 + np.sin(x1) + 1 + (1 if x0+t1*x1-t0>0 else 0)
>>> def lim0(x1, x2, t0, t1):
...     return [t0*x1 + t1*x2 - 1, t0*x1 + t1*x2 + 1]
>>> def lim1(x2, t0, t1):
...     return [x2 + t0**2*t1**3 - 1, x2 + t0**2*t1**3 + 1]
>>> def lim2(t0, t1):
...     return [t0 + t1 - 1, t0 + t1 + 1]
>>> def opts0(x1, x2, t0, t1):
...     return {'points' : [t0 - t1*x1]}
>>> def opts1(x2, t0, t1):
...     return {}
>>> def opts2(t0, t1):
...     return {}
>>> integrate.nquad(func2, [lim0, lim1, lim2], args=(0,1),
...                 opts=[opts0, opts1, opts2])
(36.099919226771625, 1.8546948553373528e-07) 

scipy.integrate.fixed_quad

原文链接:docs.scipy.org/doc/scipy-1.12.0/reference/generated/scipy.integrate.fixed_quad.html#scipy.integrate.fixed_quad

scipy.integrate.fixed_quad(func, a, b, args=(), n=5)

使用固定阶高斯积分法计算定积分。

使用高斯积分法从ab积分func,积分阶数为n

参数:

func可调用

一个 Python 函数或方法用于积分(必须接受向量输入)。如果积分的是一个向量值函数,则返回的数组必须具有形状(..., len(x))

a浮点数

积分的下限

b浮点数

积分的上限

args元组,可选

传递给函数的额外参数(如果有)。

n整数,可选

积分的阶数。默认值为 5。

返回:

val浮点数

高斯积分法对积分的近似

none

静态返回的空值

另请参阅

quad

使用 QUADPACK 的自适应积分

dblquad

双重积分

tplquad

三重积分

romberg

自适应 Romberg 积分

quadrature

自适应高斯积分

romb

采样数据的积分器

simpson

采样数据的积分器

cumulative_trapezoid

采样数据的累积积分

ode

ODE(常微分方程)积分器

odeint

ODE(常微分方程)积分器

示例

>>> from scipy import integrate
>>> import numpy as np
>>> f = lambda x: x**8
>>> integrate.fixed_quad(f, 0.0, 1.0, n=4)
(0.1110884353741496, None)
>>> integrate.fixed_quad(f, 0.0, 1.0, n=5)
(0.11111111111111102, None)
>>> print(1/9.0)  # analytical result
0.1111111111111111 
>>> integrate.fixed_quad(np.cos, 0.0, np.pi/2, n=4)
(0.9999999771971152, None)
>>> integrate.fixed_quad(np.cos, 0.0, np.pi/2, n=5)
(1.000000000039565, None)
>>> np.sin(np.pi/2)-np.sin(0)  # analytical result
1.0 

scipy.integrate.quadrature

原文:docs.scipy.org/doc/scipy-1.12.0/reference/generated/scipy.integrate.quadrature.html#scipy.integrate.quadrature

scipy.integrate.quadrature(func, a, b, args=(), tol=1.49e-08, rtol=1.49e-08, maxiter=50, vec_func=True, miniter=1)

使用固定容差高斯积分计算定积分。

自 SciPy 1.12.0 版本起已弃用:此函数已自 SciPy 1.12.0 版本起弃用,并将在 SciPy 1.15.0 版本中移除。请改用scipy.integrate.quad函数。

使用绝对容差tolab积分func的累积高斯积分。

参数:

func函数。

用于积分的 Python 函数或方法。

afloat。

积分下限。

bfloat。

积分上限。

argstuple,可选。

传递给函数的额外参数。

tol, rtolfloat,可选。

当最后两次迭代之间的误差小于tol或相对变化小于rtol时停止迭代。

maxiterint,可选。

高斯积分的最大阶数。

vec_funcbool,可选。

True 或 False 表示 func 是否处理数组作为参数(是“向量”函数)。默认为 True。

miniterint,可选。

高斯积分的最小阶数。

返回:

valfloat。

高斯积分的近似(在容差范围内)到积分。

errfloat。

积分估计的最后两次差异。

另请参阅:

romberg函数。

自适应的 Romberg 积分。

fixed_quad函数。

固定阶数的高斯积分。

quad函数。

使用 QUADPACK 进行自适应积分。

dblquad函数。

双重积分。

tplquad函数。

三重积分。

romb函数。

用于采样数据的积分器。

simpson函数。

用于采样数据的积分器。

cumulative_trapezoid函数。

用于采样数据的累积积分。

ode函数。

ODE 积分器。

odeint函数。

ODE 积分器。

示例。

>>> from scipy import integrate
>>> import numpy as np
>>> f = lambda x: x**8
>>> integrate.quadrature(f, 0.0, 1.0)
(0.11111111111111106, 4.163336342344337e-17)
>>> print(1/9.0)  # analytical result
0.1111111111111111 
>>> integrate.quadrature(np.cos, 0.0, np.pi/2)
(0.9999999999999536, 3.9611425250996035e-11)
>>> np.sin(np.pi/2)-np.sin(0)  # analytical result
1.0 

scipy.integrate.romberg

原文链接:docs.scipy.org/doc/scipy-1.12.0/reference/generated/scipy.integrate.romberg.html#scipy.integrate.romberg

scipy.integrate.romberg(function, a, b, args=(), tol=1.48e-08, rtol=1.48e-08, show=False, divmax=10, vec_func=False)

一个可调用函数或方法的 Romberg 积分。

自 SciPy 1.12.0 版开始弃用:该函数在 SciPy 1.12.0 版弃用,将在 SciPy 1.15.0 版中删除。请使用 scipy.integrate.quad 替代。

返回函数 function(一个一维变量的函数)在区间(ab)上的积分。

如果 show 设为 1,则会打印出中间结果的三角形数组。如果 vec_func 为真(默认为假),则假定 function 支持向量参数。

参数:

functioncallable

要积分的函数。

afloat

积分的下限。

bfloat

积分的上限。

返回:

resultsfloat

积分结果。

其他参数:

argstuple, 可选

要传递给函数的额外参数。每个 args 的元素将作为单个参数传递给 func。默认不传递任何额外参数。

tol, rtolfloat, 可选

所需的绝对和相对容差。默认值为 1.48e-8。

showbool, 可选

是否打印结果。默认为假。

divmaxint, 可选

最大外推阶数。默认为 10。

vec_funcbool, 可选

func 是否处理数组作为参数(即是否为“向量”函数)。默认为假。

另请参见

fixed_quad

固定阶高斯积分。

quad

使用 QUADPACK 的自适应积分。

dblquad

双重积分。

tplquad

三重积分。

romb

采样数据的积分器。

simpson

采样数据的积分器。

cumulative_trapezoid

采样数据的累积积分。

ode

ODE 积分器。

odeint

ODE 积分器。

参考文献

[1]

‘Romberg 方法’ en.wikipedia.org/wiki/Romberg%27s_method

示例

积分高斯函数从 0 到 1 并与误差函数进行比较。

>>> from scipy import integrate
>>> from scipy.special import erf
>>> import numpy as np
>>> gaussian = lambda x: 1/np.sqrt(np.pi) * np.exp(-x**2)
>>> result = integrate.romberg(gaussian, 0, 1, show=True)
Romberg integration of <function vfunc at ...> from [0, 1] 
Steps  StepSize  Results
    1  1.000000  0.385872
    2  0.500000  0.412631  0.421551
    4  0.250000  0.419184  0.421368  0.421356
    8  0.125000  0.420810  0.421352  0.421350  0.421350
   16  0.062500  0.421215  0.421350  0.421350  0.421350  0.421350
   32  0.031250  0.421317  0.421350  0.421350  0.421350  0.421350  0.421350 

最终结果是在 33 个函数评估后为 0.421350396475。

>>> print("%g  %g" % (2*result, erf(1)))
0.842701 0.842701 

scipy.integrate.newton_cotes

原文链接:docs.scipy.org/doc/scipy-1.12.0/reference/generated/scipy.integrate.newton_cotes.html#scipy.integrate.newton_cotes

scipy.integrate.newton_cotes(rn, equal=0)

返回牛顿-科特斯积分的权重和误差系数。

假设我们在位置为 x_0, x_1, …, x_N 的(N+1)个样本上有 f 的样本。那么在 x_0 和 x_N 之间的 N 点牛顿-科特斯公式为:

(\int_{x_0}^{x_N} f(x)dx = \Delta x \sum_{i=0}^{N} a_i f(x_i) + B_N (\Delta x)^{N+2} f^{N+1} (\xi))

其中 (\xi \in [x_0,x_N]),(\Delta x = \frac{x_N-x_0}{N}) 是平均样本间距。

如果样本等间隔且 N 为偶数,则误差项为 (B_N (\Delta x)^{N+3} f^{N+2}(\xi))。

参数:

rnint

整数阶等间隔数据或样本相对位置,其中第一个样本为 0,最后一个为 N,其中 N+1 为rn的长度。N 为牛顿-科特斯积分的阶数。

equalint, 可选

设为 1 以强制等间隔数据。

返回:

anndarray

1-D 权重数组,应用于提供的样本位置处的函数。

Bfloat

错误系数。

注意事项

通常,牛顿-科特斯规则用于较小的积分区域,并且使用复合规则返回总积分。

示例

计算在[0, (\pi)]内 sin(x)的积分:

>>> from scipy.integrate import newton_cotes
>>> import numpy as np
>>> def f(x):
...     return np.sin(x)
>>> a = 0
>>> b = np.pi
>>> exact = 2
>>> for N in [2, 4, 6, 8, 10]:
...     x = np.linspace(a, b, N + 1)
...     an, B = newton_cotes(N, 1)
...     dx = (b - a) / N
...     quad = dx * np.sum(an * f(x))
...     error = abs(quad - exact)
...     print('{:2d}  {:10.9f}  {:.5e}'.format(N, quad, error))
...
 2   2.094395102   9.43951e-02
 4   1.998570732   1.42927e-03
 6   2.000017814   1.78136e-05
 8   1.999999835   1.64725e-07
10   2.000000001   1.14677e-09 

scipy.integrate.qmc_quad

原文链接:docs.scipy.org/doc/scipy-1.12.0/reference/generated/scipy.integrate.qmc_quad.html#scipy.integrate.qmc_quad

scipy.integrate.qmc_quad(func, a, b, *, n_estimates=8, n_points=1024, qrng=None, log=False)

使用准蒙特卡洛积分计算 N 维积分。

参数:

func可调用对象

积分被积函数。必须接受单个参数x,一个数组,指定要评估标量值积分被积函数的点,并返回被积函数的值。为了效率,该函数应该向量化,接受形状为(d, n_points)的数组,其中d是变量的数量(即函数域的维度),n_points是积分点的数量,返回形状为(n_points,)的数组,即每个积分点的被积函数值。

a, b类数组

一维数组,分别指定d个变量的积分下限和上限。

n_estimates, n_points整数,可选

n_estimates(默认值:8)统计独立的 QMC 样本,每个n_points(默认值:1024)点,将由qrng生成。函数func将在n_points * n_estimates个点上进行评估。详见注释。

qrngQMCEngine,可选

QMCEngine 的实例,用于抽样 QMC 点。QMCEngine 必须初始化为与传递给func的变量x1, ..., xd对应的维数d。提供的 QMCEngine 用于生成第一个积分估计值。如果n_estimates大于 1,则从第一个 QMCEngine 生成额外的 QMCEngine(如果有选项则启用混淆)。如果未提供 QMCEngine,则将使用默认的scipy.stats.qmc.Halton,其维数由a的长度确定。

log布尔值,默认值:False

当设置为 True 时,func返回积分被积函数的对数,结果对象包含积分的对数。

返回:

result对象

具有以下属性的结果对象:

积分值浮点数

积分估计值。

standard_error:

误差估计。详见注释以获取解释。

注释

在 QMC 样本的 n_points 点中的积分值被用来产生对积分的估计。这个估计来自于可能的积分估计的一个群体,我们获得的值取决于评估积分的特定点。我们对此过程执行 n_estimates 次,每次评估不同混乱的 QMC 点的积分值,有效地从积分估计的群体中抽取 i.i.d. 随机样本。这些积分估计的样本均值 (m) 是真实积分值的无偏估计,而这些估计的样本均值 (s) 的标准误差可以使用自由度为 n_estimates - 1 的 t 分布生成置信区间。或许反直觉地,增加 n_points 而保持总的函数评估点数 n_points * n_estimates 固定倾向于减少实际误差,而增加 n_estimates 则倾向于减少误差估计。

示例

QMC(低差异序列蒙特卡罗)求积法在计算高维积分时特别有用。一个例子积分被用作多元正态分布的概率密度函数。

>>> import numpy as np
>>> from scipy import stats
>>> dim = 8
>>> mean = np.zeros(dim)
>>> cov = np.eye(dim)
>>> def func(x):
...     # `multivariate_normal` expects the _last_ axis to correspond with
...     # the dimensionality of the space, so `x` must be transposed
...     return stats.multivariate_normal.pdf(x.T, mean, cov) 

要计算单位超立方体上的积分:

>>> from scipy.integrate import qmc_quad
>>> a = np.zeros(dim)
>>> b = np.ones(dim)
>>> rng = np.random.default_rng()
>>> qrng = stats.qmc.Halton(d=dim, seed=rng)
>>> n_estimates = 8
>>> res = qmc_quad(func, a, b, n_estimates=n_estimates, qrng=qrng)
>>> res.integral, res.standard_error
(0.00018429555666024108, 1.0389431116001344e-07) 

对积分的双边、99% 置信区间可以估计为:

>>> t = stats.t(df=n_estimates-1, loc=res.integral,
...             scale=res.standard_error)
>>> t.interval(0.99)
(0.0001839319802536469, 0.00018465913306683527) 

确实,scipy.stats.multivariate_normal 返回的数值在这个范围内。

>>> stats.multivariate_normal.cdf(b, mean, cov, lower_limit=a)
0.00018430867675187443 

scipy.integrate.IntegrationWarning

原文:docs.scipy.org/doc/scipy-1.12.0/reference/generated/scipy.integrate.IntegrationWarning.html#scipy.integrate.IntegrationWarning

exception scipy.integrate.IntegrationWarning

整合过程中的问题警告。

with_traceback()

Exception.with_traceback(tb) – 设置 self.traceback 为 tb 并返回 self。

scipy.integrate.AccuracyWarning

原文:docs.scipy.org/doc/scipy-1.12.0/reference/generated/scipy.integrate.AccuracyWarning.html#scipy.integrate.AccuracyWarning

exception scipy.integrate.AccuracyWarning
with_traceback()

Exception.with_traceback(tb) – 设置 self.traceback 为 tb 并返回 self。

scipy.integrate.trapezoid

原文:docs.scipy.org/doc/scipy-1.12.0/reference/generated/scipy.integrate.trapezoid.html#scipy.integrate.trapezoid

scipy.integrate.trapezoid(y, x=None, dx=1.0, axis=-1)

使用复合梯形法则沿给定轴积分。

如果提供了x,则积分按其元素的顺序进行 - 它们未排序。

沿给定轴上的每个 1d 切片积分y(x),计算(\int y(x) dx)。当指定x时,这将沿参数曲线积分,计算(\int_t y(t) dt = \int_t y(t) \left.\frac{dx}{dt}\right|_{x=x(t)} dt)。

参数:

yarray_like

输入要积分的数组。

xarray_like,可选

对应于y值的样本点。如果x为 None,则假定样本点等间隔地间隔dx。默认为 None。

dx标量,可选

x为 None 时,样本点之间的间距。默认为 1。

axisint,可选

要积分的轴。

返回:

trapezoidfloat 或 ndarray

定义y = n 维数组的定积分,通过梯形法则沿单轴近似。如果y是一维数组,则结果为浮点数。如果n大于 1,则结果是一个n-1 维数组。

另见

cumulative_trapezoid, simpson, romb

注意事项

图像[2]说明梯形法则 - y 轴点的位置将从y数组中获取,默认情况下,点之间的 x 轴距离将为 1.0,也可以使用x数组或dx标量提供。返回值将等于红线下的联合面积。

参考

[1]

Wikipedia 页面:en.wikipedia.org/wiki/Trapezoidal_rule

[2]

插图:en.wikipedia.org/wiki/File:Composite_trapezoidal_rule_illustration.png

示例

应用梯形法则在均匀间隔的点上:

>>> import numpy as np
>>> from scipy import integrate
>>> integrate.trapezoid([1, 2, 3])
4.0 

可以通过xdx参数选择样本点之间的间距:

>>> integrate.trapezoid([1, 2, 3], x=[4, 6, 8])
8.0
>>> integrate.trapezoid([1, 2, 3], dx=2)
8.0 

使用递减的x相当于反向积分:

>>> integrate.trapezoid([1, 2, 3], x=[8, 6, 4])
-8.0 

更一般地,使用x来沿参数曲线积分。我们可以估计积分(\int_0¹ x² = 1/3):

>>> x = np.linspace(0, 1, num=50)
>>> y = x**2
>>> integrate.trapezoid(y, x)
0.33340274885464394 

或者估计圆的面积,注意我们重复了闭合曲线的样本:

>>> theta = np.linspace(0, 2 * np.pi, num=1000, endpoint=True)
>>> integrate.trapezoid(np.cos(theta), x=np.sin(theta))
3.141571941375841 

trapezoid可以沿指定轴应用以在一次调用中进行多个计算:

>>> a = np.arange(6).reshape(2, 3)
>>> a
array([[0, 1, 2],
 [3, 4, 5]])
>>> integrate.trapezoid(a, axis=0)
array([1.5, 2.5, 3.5])
>>> integrate.trapezoid(a, axis=1)
array([2.,  8.]) 

scipy.integrate.cumulative_trapezoid

原文链接:docs.scipy.org/doc/scipy-1.12.0/reference/generated/scipy.integrate.cumulative_trapezoid.html#scipy.integrate.cumulative_trapezoid

scipy.integrate.cumulative_trapezoid(y, x=None, dx=1.0, axis=-1, initial=None)

使用复合梯形法累积地对 y(x) 进行积分。

参数:

y array_like

要积分的值。

x array_like,可选

要进行积分的坐标。如果为 None(默认),则在 y 的连续元素之间使用 dx 的间距。

dx 浮点数,可选

y 元素之间的间距。仅当 x 为 None 时使用。

axis 整数,可选

指定要累积的轴。默认为 -1(最后一个轴)。

initial 标量,可选

如果给定,将此值插入返回结果的开头。仅接受 0 或 None。默认为 None,这意味着 res 沿积分轴比 y 少一个元素。

自 SciPy 1.15.0 版本起已弃用:initial 的非零输入选项。此后,如果 initial 不为 None 或 0,则会引发 ValueError。

返回:

res ndarray

y 沿 axis 累积积分的结果。如果 initial 为 None,则形状使得积分轴比 y 的轴少一个值。如果给定了 initial,则形状与 y 相同。

另请参阅

numpy.cumsum, numpy.cumprod

cumulative_simpson

使用 Simpson's 1/3 规则的累积积分

quad

使用 QUADPACK 的自适应积分

romberg

自适应 Romberg 积分

quadrature

自适应高斯积分

fixed_quad

固定阶数的高斯积分

dblquad

双重积分

tplquad

三重积分

romb

用于采样数据的积分器

ode

ODE 积分器

odeint

ODE 积分器

示例

>>> from scipy import integrate
>>> import numpy as np
>>> import matplotlib.pyplot as plt 
>>> x = np.linspace(-2, 2, num=20)
>>> y = x
>>> y_int = integrate.cumulative_trapezoid(y, x, initial=0)
>>> plt.plot(x, y_int, 'ro', x, y[0] + 0.5 * x**2, 'b-')
>>> plt.show() 

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

scipy.integrate.simpson

原文:docs.scipy.org/doc/scipy-1.12.0/reference/generated/scipy.integrate.simpson.html#scipy.integrate.simpson

scipy.integrate.simpson(y, *, x=None, dx=1.0, axis=-1, even=<object object>)

使用给定轴上的样本和复合 Simpson 规则来积分 y(x)。如果 x 为 None,则假定 dx 的间距。

如果有偶数个样本 N,则有奇数个间隔(N-1),但 Simpson 规则需要偶数个间隔。参数‘even’控制如何处理此问题。

参数:

yarray_like

被积数组。

xarray_like,可选

如果给定,则为y进行采样的点。

dxfloat,可选

沿x轴的积分点间距。仅当x为 None 时使用。默认为 1。

axisint,可选

进行积分的轴。默认为最后一个轴。

even,可选

‘avg’平均两个结果:

  1. 使用第一个 N-2 个间隔和最后一个间隔上的梯形法则。

  2. 使用最后 N-2 个间隔和第一个间隔上的梯形法则。

‘first’对前 N-2 个间隔使用 Simpson 规则

最后一个间隔上的梯形法则。

‘last’对最后 N-2 个间隔使用 Simpson 规则进行

第一个间隔上的梯形法则。

None:等同于‘simpson’(默认)

‘simpson’使用 Simpson 规则对前 N-2 个间隔进行积分。

添加一个由 Cartwright[1]提出的 3 点抛物线段到最后一个间隔中。如果要积分的轴只有两个点,则积分回退到梯形积分。

版本 1.11.0 中的新功能。

从版本 1.11.0 开始更改:新添加的‘simpson’选项现在是默认选项,因为在大多数情况下更准确。

自版本 1.11.0 起弃用:参数even已弃用,并将在 SciPy 1.14.0 中删除。此后,偶数点数的行为将遵循even='simpson'

返回:

浮点数。

使用复合 Simpson 规则计算的估计积分。

另请参见

quad

使用 QUADPACK 进行自适应积分。

romberg

自适应 Romberg 积分。

quadrature

自适应高斯积分。

fixed_quad

固定顺序的高斯积分。

dblquad

双重积分。

tplquad

三重积分。

romb

用于采样数据的积分器

cumulative_trapezoid

用于采样数据的累积积分

cumulative_simpson

使用 Simpson’s 1/3 规则进行累积积分

ode

ODE(常微分方程)积分器

odeint

ODE(常微分方程)积分器

笔记

对于等间隔的样本数目为奇数的情况,如果函数是三阶或更低阶的多项式,则结果是精确的。如果样本不是等间隔的,则结果仅在函数为二阶或更低阶的多项式时是精确的。

参考文献

[1]

Kenneth V. Cartwright. 使用 MS Excel 和不规则间隔数据的 Simpson’s Rule Cumulative Integration。《数学科学与数学教育杂志》。12 (2): 1-9

示例

>>> from scipy import integrate
>>> import numpy as np
>>> x = np.arange(0, 10)
>>> y = np.arange(0, 10) 
>>> integrate.simpson(y, x)
40.5 
>>> y = np.power(x, 3)
>>> integrate.simpson(y, x)
1640.5
>>> integrate.quad(lambda x: x**3, 0, 9)[0]
1640.25 
>>> integrate.simpson(y, x, even='first')
1644.5 

scipy.integrate.cumulative_simpson

原文链接:docs.scipy.org/doc/scipy-1.12.0/reference/generated/scipy.integrate.cumulative_simpson.html#scipy.integrate.cumulative_simpson

scipy.integrate.cumulative_simpson(y, *, x=None, dx=1.0, axis=-1, initial=None)

使用复合辛普森 1/3 法累积积分y(x)。假定每个点与其两个相邻点之间存在二次关系来计算每个点的积分样本。

参数:

yarray_like

需要积分的值。至少需要沿着的一个点。如果提供的点少于或等于两个,则不可能使用辛普森积分法,结果将使用cumulative_trapezoid计算。

x数组型,可选

要进行积分的坐标。必须与y具有相同形状或在上具有与y相同长度的 1D 数组。x还必须在上严格递增。如果x为 None(默认),则使用y中连续元素之间的间距dx进行积分。

dx标量或数组型,可选

y元素之间的间距。仅在x为 None 时使用。可以是浮点数,也可以是与y相同形状但在上长度为一的数组。默认为 1.0。

axis整数,可选

指定要沿其进行积分的。默认为-1(最后一个轴)。

initial标量或数组型,可选

如果提供,则在返回结果的开头插入该值,并将其添加到其余结果中。默认为 None,这意味着在x[0]处不返回任何值,并且res沿积分轴比y少一个元素。可以是浮点数,也可以是与y相同形状但在上长度为一的数组。

返回:

resndarray

沿积分y的累积结果。如果initial为 None,则形状使得积分轴比y少一个值。如果给定initial,形状与y相同。

另请参阅

numpy.cumsum

cumulative_trapezoid

使用复合梯形法进行累积积分

simpson

用于采样数据的复合辛普森法积分器

自 1.12.0 版开始新引入。

复合辛普森 1/3 法可用于近似采样输入函数y(x)的定积分 [1]。该方法假定在包含任意三个连续采样点的区间上存在二次关系。

考虑三个连续点:((x_1, y_1), (x_2, y_2), (x_3, y_3))。

假设在这三个点上存在二次关系,(x_1)和(x_2)之间的子区间积分由[2]的公式(8)给出:

[\begin{split}\int_{x_1}^{x_2} y(x) dx\ &= \frac{x_2-x_1}{6}\left[\ \left{3-\frac{x_2-x_1}{x_3-x_1}\right} y_1 + \ \left{3 + \frac{(x_2-x_1)²}{(x_3-x_2)(x_3-x_1)} + \ \frac{x_2-x_1}{x_3-x_1}\right} y_2\ - \frac{(x_2-x_1)²}{(x_3-x_2)(x_3-x_1)} y_3\right]\end{split}]

在(x_2)和(x_3)之间的积分通过交换(x_1)和(x_3)的位置来计算。对每个子区间分别进行估计积分,然后累加以获得最终结果。

对于等间距样本,如果函数是三次或更低次的多项式,并且子区间数是偶数,则结果是精确的[1]。否则,对于二次或更低次的多项式,积分是精确的。

参考文献

[1] (1,2)

Wikipedia 页面:en.wikipedia.org/wiki/Simpson’s_rule

[2]

Cartwright, Kenneth V. Simpson’s Rule Cumulative Integration with MS Excel and Irregularly-spaced Data. Journal of Mathematical Sciences and Mathematics Education. 12 (2): 1-9

示例

>>> from scipy import integrate
>>> import numpy as np
>>> import matplotlib.pyplot as plt
>>> x = np.linspace(-2, 2, num=20)
>>> y = x**2
>>> y_int = integrate.cumulative_simpson(y, x=x, initial=0)
>>> fig, ax = plt.subplots()
>>> ax.plot(x, y_int, 'ro', x, x**3/3 - (x[0])**3/3, 'b-')
>>> ax.grid()
>>> plt.show() 

../../_images/scipy-integrate-cumulative_simpson-1_00_00.png

cumulative_simpson 的输出类似于连续调用 simpson,每次的积分上限逐渐增加,但并非完全相同。

>>> def cumulative_simpson_reference(y, x):
...     return np.asarray([integrate.simpson(y[:i], x=x[:i])
...                        for i in range(2, len(y) + 1)])
>>>
>>> rng = np.random.default_rng()
>>> x, y = rng.random(size=(2, 10))
>>> x.sort()
>>>
>>> res = integrate.cumulative_simpson(y, x=x)
>>> ref = cumulative_simpson_reference(y, x)
>>> equal = np.abs(res - ref) < 1e-15
>>> equal  # not equal when `simpson` has even number of subintervals
array([False,  True, False,  True, False,  True, False,  True,  True]) 

这是预期的结果:因为 cumulative_simpson 拥有比 simpson 更多的信息,通常可以在子区间上产生更精确的积分估计。

scipy.integrate.romb

原文链接:docs.scipy.org/doc/scipy-1.12.0/reference/generated/scipy.integrate.romb.html#scipy.integrate.romb

scipy.integrate.romb(y, dx=1.0, axis=-1, show=False)

使用函数的样本进行 Romberg 积分。

参数:

yarray_like

函数的2**k + 1等间距样本的向量。

dxfloat,可选

样本间距。默认为 1。

axisint,可选

要进行积分的轴。默认为-1(最后一个轴)。

showbool,可选

y是单个 1-D 数组时,如果此参数为 True,则打印从样本中的 Richardson 外推的表格。默认为 False。

返回:

rombndarray

axis的整合结果。

另请参阅

quad

使用 QUADPACK 的自适应积分

romberg

自适应的 Romberg 积分

quadrature

自适应的高斯积分

fixed_quad

固定顺序的高斯积分

dblquad

双重积分

tplquad

三重积分

simpson

用于采样数据的积分器

cumulative_trapezoid

对采样数据的累积积分

ode

ODE(常微分方程)积分器

odeint

ODE(常微分方程)积分器

示例

>>> from scipy import integrate
>>> import numpy as np
>>> x = np.arange(10, 14.25, 0.25)
>>> y = np.arange(3, 12) 
>>> integrate.romb(y)
56.0 
>>> y = np.sin(np.power(x, 2.5))
>>> integrate.romb(y)
-0.742561336672229 
>>> integrate.romb(y, show=True)
Richardson Extrapolation Table for Romberg Integration
======================================================
-0.81576
 4.63862  6.45674
-1.10581 -3.02062 -3.65245
-2.57379 -3.06311 -3.06595 -3.05664
-1.34093 -0.92997 -0.78776 -0.75160 -0.74256
======================================================
-0.742561336672229  # may vary 

scipy.integrate.solve_ivp

原文链接:docs.scipy.org/doc/scipy-1.12.0/reference/generated/scipy.integrate.solve_ivp.html#scipy.integrate.solve_ivp

scipy.integrate.solve_ivp(fun, t_span, y0, method='RK45', t_eval=None, dense_output=False, events=None, vectorized=False, args=None, **options)

为一组 ODE 解决初始值问题。

此函数对给定初始值数值积分一组普通微分方程系统:

dy / dt = f(t, y)
y(t0) = y0 

这里 t 是一维独立变量(时间),y(t)是一个 N 维向量值函数(状态),而一个 N 维向量值函数 f(t, y)确定了微分方程。目标是找到一个近似满足微分方程的 y(t),给定初始值 y(t0)=y0。

一些求解器支持在复数域中进行积分,但请注意,对于刚性 ODE 求解器,右手边必须是复可微的(满足柯西-黎曼方程[11])。要解决复数域中的问题,请使用具有复数数据类型的y0。另一个始终可用的选项是将问题分别重写为实部和虚部。

参数:

fun可调用函数

系统的右手边:在时间 t 处状态y的时间导数。调用签名是fun(t, y),其中t是标量,y是一个长度为len(y0)的 ndarray。如果使用了args(参见args参数的文档),需要传递额外的参数。fun必须返回与y相同形状的数组。详见向量化以获取更多信息。

t_span的 2 元序列

积分区间(t0, tf)。求解器从 t=t0 开始积分,直到达到 t=tf。t0 和 tf 都必须是浮点数或可以转换为浮点数的值。

y0array_like,形状(n,)

初始状态。对于复数域中的问题,请使用具有复数数据类型的y0(即使初始值是纯实数)。

method字符串或OdeSolver,可选

要使用的积分方法:

  • ‘RK45’(默认):五阶(四阶)显式 Runge-Kutta 方法[1]。通过假定四阶方法的准确性来控制误差,但使用五阶准确公式(进行本地外推)。对于密集输出,使用四次插值多项式[2]。可应用于复数域。
  • ‘RK23’:三阶(二阶)显式 Runge-Kutta 方法[3]。通过假定二阶方法的准确性来控制误差,但使用三阶准确公式(进行本地外推)。对于密集输出,使用三次 Hermite 多项式。可应用于复数域。
  • ‘DOP853’:8 阶显式 Runge-Kutta 方法[13]。Python 实现了最初在 Fortran 中编写的“DOP853”算法[14]。用于密集输出的 7 阶插值多项式精确到 7 阶。可应用于复数域。
  • ‘Radau’:Radau IIA 族的隐式 Runge-Kutta 方法,阶数为 5[4]。误差通过三阶准确的嵌入公式控制。满足配点条件的三次多项式用于密集输出。
  • ‘BDF’:基于后向差分公式的隐式多步变阶(1 到 5 阶)方法,用于导数近似[5]。该实现遵循[6]中描述的方法。采用准恒定步长方案,并使用 NDF 修正增强精度。可应用于复数域。
  • ‘LSODA’:Adams/BDF 方法,具有自动刚性检测和切换[7][8]。这是 ODEPACK 中 Fortran 求解器的包装器。

显式 Runge-Kutta 方法(‘RK23’、‘RK45’、‘DOP853’)应用于非刚性问题,而隐式方法(‘Radau’、‘BDF’)应用于刚性问题[9]。在 Runge-Kutta 方法中,推荐使用‘DOP853’以实现高精度求解(低rtolatol值)。

如果不确定,请先尝试运行‘RK45’。如果它执行异常多的迭代、发散或失败,则您的问题可能是刚性的,您应该使用‘Radau’或‘BDF’。‘LSODA’也可以是一个好的通用选择,但它可能不太方便使用,因为它包装了旧的 Fortran 代码。

您还可以传递从OdeSolver派生的任意类,该类实现了求解器。

t_eval array_like 或 None,可选

计算解的时间点,必须按顺序排序并位于t_span内。如果为 None(默认),则使用求解器选择的点。

dense_output bool,可选

是否计算连续解。默认为 False。

events callable 或 callable 列表,可选

跟踪事件。如果为 None(默认),将不会跟踪任何事件。每个事件发生在时间和状态的连续函数的零点上。每个函数必须具有event(t, y)的签名,在使用args(参见args参数的文档)时必须传递额外的参数。每个函数必须返回一个浮点数。求解器将使用根查找算法找到使得event(t, y(t)) = 0的准确时间t。默认情况下,将找到所有的零点。求解器在每一步中寻找符号变化,因此如果在一步内发生多次零点穿越,则可能会错过事件。此外,每个event函数可能具有以下属性:

terminal:bool,可选

是否在此事件发生时终止积分。如果未分配,则隐式为 False。

direction: float, optional

零点穿越的方向。如果direction为正,当从负到正时event将触发,如果direction为负,则反之。如果为 0,则任何方向都会触发事件。如果未指定,则隐式为 0。

您可以为 Python 中的任何函数分配属性,例如event.terminal = True

vectorizedbool, optional

是否可以以向量化方式调用fun。默认为 False。

如果vectorized为 False,则fun将始终使用形状为(n,)y调用,其中n = len(y0)

如果vectorized为 True,则可以使用形状为(n, k)y调用fun,其中k是整数。在这种情况下,fun必须表现出fun(t, y)[:, i] == fun(t, y[:, i])(即返回数组的每一列都是与y的每一列对应的状态的时间导数)。

设置vectorized=True允许通过‘Radau’和‘BDF’方法更快地进行雅可比矩阵的有限差分逼近,但会导致其他方法和在某些情况下(例如小的len(y0))‘Radau’和‘BDF’方法的执行速度变慢。

argstuple,可选

传递给用户定义函数的附加参数。如果给出,则所有用户定义函数都会传递这些附加参数。例如,如果fun的签名为fun(t, y, a, b, c),则jac(如果给出)和任何事件函数必须具有相同的签名,并且args必须是长度为 3 的元组。

**options

传递给选择的求解器的选项。列出了所有已实现求解器可用的选项。

first_stepfloat 或 None,可选

初始步长。默认为None,表示算法应选择。

max_stepfloat, optional

允许的最大步长。默认为 np.inf,即步长不受限制,仅由求解器确定。

rtol, atolfloat 或 array_like,可选

相对和绝对容差。求解器保持局部误差估计小于atol + rtol * abs(y)。其中rtol控制相对精度(正确数字的数量),而atol控制绝对精度(正确小数位数)。为了实现所需的rtol,请将atol设置为小于从rtol * abs(y)可以期望的最小值,以便rtol主导允许的误差。如果atol大于rtol * abs(y),则不能保证正确数字的数量。相反,为了实现所需的atol,请设置rtol,使得rtol * abs(y)始终小于atol。如果 y 的组件具有不同的比例,则通过传递形状为(n,)的 array_like 为atol的不同组件设置不同的值可能是有益的。默认值为rtol为 1e-3 和atol为 1e-6。

jacarray_like、稀疏矩阵、可调用对象或 None,可选

y 系统右手边的雅可比矩阵,要求 ‘Radau’、‘BDF’ 和 ‘LSODA’ 方法。雅可比矩阵的形状为 (n, n),其元素 (i, j) 等于 d f_i / d y_j。有三种方法来定义雅可比矩阵:

  • 如果是 array_likesparse_matrix,则假定雅可比矩阵是常数。不支持‘LSODA’。
  • 如果是可调用的,则假定雅可比矩阵依赖于 ty;将按需调用为 jac(t, y)。如果使用了 args(请参阅 args 参数的文档),还必须传递额外参数。对于‘Radau’和‘BDF’方法,返回值可能是稀疏矩阵。
  • 如果为 None(默认),雅可比矩阵将通过有限差分逼近。

通常建议提供雅可比矩阵而不是依赖有限差分逼近。

jac_sparsityarray_like,稀疏矩阵或 None,可选

定义雅可比矩阵的稀疏结构,用于有限差分逼近。其形状必须为 (n, n)。如果 jac 不为 None,则忽略此参数。如果雅可比矩阵每行只有几个非零元素,提供稀疏结构将极大加速计算 [10]。零条目意味着雅可比矩阵中相应的元素始终为零。如果为 None(默认),则假定雅可比矩阵是密集的。‘LSODA’不支持,见 lbanduband

lband, uband:整数或 None,可选

定义“LSODA”方法雅可比矩阵带宽的参数,即,jac[i, j] != 0 仅当 i - lband <= j <= i + uband。默认为 None。设置这些参数需要你的雅可比例程以打包格式返回雅可比矩阵:返回的数组必须有 n 列和 uband + lband + 1 行,其中雅可比矩阵对角线写入。具体而言,jac_packed[uband + i - j, j] = jac[i, j]scipy.linalg.solve_banded 中也使用相同格式(查看示例)。这些参数也可以与 jac=None 一起使用,以减少通过有限差分估计的雅可比元素数量。

min_step:浮点数,可选

“LSODA”方法的最小允许步长。默认情况下 min_step 为零。

返回:

包对象,具有以下字段定义:

tndarray,形状为 (n_points,)

时间点。

yndarray,形状为 (n, n_points)

t 处解的值。

solOdeSolutionNone

找到的解作为 OdeSolution 实例;如果 dense_output 设置为 False,则为 None

t_eventsndarray 列表或 None

包含每种事件类型检测到事件的数组列表。如果 eventsNone,则为 None

y_eventsndarray 列表或 None

对于每个t_events的值,对应的解的值。如果events为 None,则为 None。

nfev整数

右手边求值的数量。

njev整数

雅可比矩阵的求值数量。

nlu整数

LU 分解的数量。

状态整数

算法终止的原因:

  • -1:积分步骤失败。
  • 0:求解器成功到达tspan的末尾。
  • 1:发生终止事件。

消息字符串

算法终止原因的可读描述。

成功布尔值

如果求解器达到时间间隔的结束或发生终止事件(status >= 0),则为 True。

参考文献

[1]

J. R. Dormand, P. J. Prince, “一族嵌入 Runge-Kutta 公式”,计算和应用数学杂志,第 6 卷,第 1 期,pp. 19-26,1980 年。

[2]

L. W. Shampine, “一些实用的 Runge-Kutta 公式”,计算数学,第 46 卷,第 173 期,pp. 135-150,1986 年。

[3]

P. Bogacki, L.F. Shampine, “一对 3(2)的 Runge-Kutta 公式”,应用数学通讯,第 2 卷,第 4 期,pp. 321-325,1989 年。

[4]

E. Hairer, G. Wanner, “求解普通微分方程 II:刚性和微分代数问题”,第 IV.8 节。

[5]

反向差分公式在维基百科上。

[6]

L. F. Shampine, M. W. Reichelt, “MATLAB ODE 套件”,SIAM J. SCI. COMPUTE.,第 18 卷,第 1 期,pp. 1-22,1997 年。

[7]

A. C. Hindmarsh, “ODEPACK,一组系统化的常微分方程求解器”,IMACS 科学计算交易,第 1 卷,pp. 55-64,1983 年。

[8]

L. Petzold, “自动选择求解刚性和非刚性常微分方程方法”,SIAM 科学与统计计算期刊,第 4 卷,第 1 期,pp. 136-148,1983 年。

[9]

刚性方程在维基百科上。

[10]

A. Curtis, M. J. D. Powell, and J. Reid, “关于稀疏雅可比矩阵估计的研究”,数学应用研究所杂志,第 13 卷,pp. 117-120,1974 年。

[11]

柯西-黎曼方程在维基百科上。

[12]

Lotka-Volterra 方程组在维基百科上。

[13]

E. Hairer, S. P. Norsett G. Wanner, “求解普通微分方程 I:非刚性问题”,第 II 节。

[14]

DOP853 的原始 Fortran 代码页面

例子

显示自动选择时间点的基本指数衰减。

>>> import numpy as np
>>> from scipy.integrate import solve_ivp
>>> def exponential_decay(t, y): return -0.5 * y
>>> sol = solve_ivp(exponential_decay, [0, 10], [2, 4, 8])
>>> print(sol.t)
[ 0\.          0.11487653  1.26364188  3.06061781  4.81611105  6.57445806
 8.33328988 10\.        ]
>>> print(sol.y)
[[2\.         1.88836035 1.06327177 0.43319312 0.18017253 0.07483045
 0.03107158 0.01350781]
 [4\.         3.7767207  2.12654355 0.86638624 0.36034507 0.14966091
 0.06214316 0.02701561]
 [8\.         7.5534414  4.25308709 1.73277247 0.72069014 0.29932181
 0.12428631 0.05403123]] 

指定希望获得解的点。

>>> sol = solve_ivp(exponential_decay, [0, 10], [2, 4, 8],
...                 t_eval=[0, 1, 2, 4, 10])
>>> print(sol.t)
[ 0  1  2  4 10]
>>> print(sol.y)
[[2\.         1.21305369 0.73534021 0.27066736 0.01350938]
 [4\.         2.42610739 1.47068043 0.54133472 0.02701876]
 [8\.         4.85221478 2.94136085 1.08266944 0.05403753]] 

炮弹向上发射,撞击时有终端事件。通过猴子修补一个函数来应用事件的terminaldirection字段。这里的y[0]是位置,y[1]是速度。抛射物从位置 0 以+10 的速度开始。注意,积分从未达到 t=100,因为事件是终端的。

>>> def upward_cannon(t, y): return [y[1], -0.5]
>>> def hit_ground(t, y): return y[0]
>>> hit_ground.terminal = True
>>> hit_ground.direction = -1
>>> sol = solve_ivp(upward_cannon, [0, 100], [0, 10], events=hit_ground)
>>> print(sol.t_events)
[array([40.])]
>>> print(sol.t)
[0.00000000e+00 9.99900010e-05 1.09989001e-03 1.10988901e-02
 1.11088891e-01 1.11098890e+00 1.11099890e+01 4.00000000e+01] 

使用dense_outputevents找到炮弹轨迹顶点的位置,即 100。顶点并非终点,因此同时找到顶点和地面触地点。在 t=20 时没有信息,因此使用 sol 属性评估解。通过设置dense_output=True返回 sol 属性。另外,y_events属性可用于访问事件发生时的解。

>>> def apex(t, y): return y[1]
>>> sol = solve_ivp(upward_cannon, [0, 100], [0, 10],
...                 events=(hit_ground, apex), dense_output=True)
>>> print(sol.t_events)
[array([40.]), array([20.])]
>>> print(sol.t)
[0.00000000e+00 9.99900010e-05 1.09989001e-03 1.10988901e-02
 1.11088891e-01 1.11098890e+00 1.11099890e+01 4.00000000e+01]
>>> print(sol.sol(sol.t_events[1][0]))
[100\.   0.]
>>> print(sol.y_events)
[array([[-5.68434189e-14, -1.00000000e+01]]),
 array([[1.00000000e+02, 1.77635684e-15]])] 

作为具有额外参数系统的示例,我们将实现 Lotka-Volterra 方程[12]

>>> def lotkavolterra(t, z, a, b, c, d):
...     x, y = z
...     return [a*x - b*x*y, -c*y + d*x*y]
... 

我们通过args参数传入参数值 a=1.5, b=1, c=3 和 d=1。

>>> sol = solve_ivp(lotkavolterra, [0, 15], [10, 5], args=(1.5, 1, 3, 1),
...                 dense_output=True) 

计算密集解并绘制。

>>> t = np.linspace(0, 15, 300)
>>> z = sol.sol(t)
>>> import matplotlib.pyplot as plt
>>> plt.plot(t, z.T)
>>> plt.xlabel('t')
>>> plt.legend(['x', 'y'], shadow=True)
>>> plt.title('Lotka-Volterra System')
>>> plt.show() 

../../_images/scipy-integrate-solve_ivp-1_00_00.png

几个使用 solve_ivp 解决带有复杂矩阵A的微分方程y' = Ay的示例。

>>> A = np.array([[-0.25 + 0.14j, 0, 0.33 + 0.44j],
...               [0.25 + 0.58j, -0.2 + 0.14j, 0],
...               [0, 0.2 + 0.4j, -0.1 + 0.97j]]) 

用上述Ay作为 3x1 向量解决 IVP:

>>> def deriv_vec(t, y):
...     return A @ y
>>> result = solve_ivp(deriv_vec, [0, 25],
...                    np.array([10 + 0j, 20 + 0j, 30 + 0j]),
...                    t_eval=np.linspace(0, 25, 101))
>>> print(result.y[:, 0])
[10.+0.j 20.+0.j 30.+0.j]
>>> print(result.y[:, -1])
[18.46291039+45.25653651j 10.01569306+36.23293216j
 -4.98662741+80.07360388j] 

用上面的A解决 IVP,其中y是 3x3 矩阵:

>>> def deriv_mat(t, y):
...     return (A @ y.reshape(3, 3)).flatten()
>>> y0 = np.array([[2 + 0j, 3 + 0j, 4 + 0j],
...                [5 + 0j, 6 + 0j, 7 + 0j],
...                [9 + 0j, 34 + 0j, 78 + 0j]]) 
>>> result = solve_ivp(deriv_mat, [0, 25], y0.flatten(),
...                    t_eval=np.linspace(0, 25, 101))
>>> print(result.y[:, 0].reshape(3, 3))
[[ 2.+0.j  3.+0.j  4.+0.j]
 [ 5.+0.j  6.+0.j  7.+0.j]
 [ 9.+0.j 34.+0.j 78.+0.j]]
>>> print(result.y[:, -1].reshape(3, 3))
[[  5.67451179 +12.07938445j  17.2888073  +31.03278837j
 37.83405768 +63.25138759j]
 [  3.39949503 +11.82123994j  21.32530996 +44.88668871j
 53.17531184+103.80400411j]
 [ -2.26105874 +22.19277664j -15.1255713  +70.19616341j
 -38.34616845+153.29039931j]] 

scipy.integrate.RK23

原文:docs.scipy.org/doc/scipy-1.12.0/reference/generated/scipy.integrate.RK23.html#scipy.integrate.RK23

class scipy.integrate.RK23(fun, t0, y0, t_bound, max_step=inf, rtol=0.001, atol=1e-06, vectorized=False, first_step=None, **extraneous)

显式三阶 Runge-Kutta 方法(2 阶)。

这使用 Bogacki-Shampine 配对的公式[1]。误差受二阶方法精度控制,但使用三阶准确公式进行步骤(局部外推完成)。稠密输出使用三次 Hermite 多项式。

可应用于复数域。

参数:

fun可调用对象

系统的右手边:在时间t处状态y的时间导数。调用签名为fun(t, y),其中t是标量,y是具有len(y) = len(y0)形状的 ndarray。fun必须返回与y相同形状的数组。有关更多信息,请参见vectorized

t0浮点数

初始时间。

y0array_like,形状为(n,)

初始状态。

t_bound浮点数

边界时间 - 积分不会超出此时间。它还确定积分的方向。

first_step浮点数或 None,可选

初始步长。默认为None,表示算法应选择。

max_step浮点数,可选

允许的最大步长。默认为 np.inf,即步长无界,完全由求解器确定。

rtol, atol浮点数和 array_like,可选

相对和绝对容差。求解器保持局部误差估计小于atol + rtol * abs(y)。这里rtol控制相对精度(正确数字的数量),而atol控制绝对精度(正确小数位数)。为了达到期望的rtol,将atol设置为比rtol * abs(y)预期的最小值更小,以便rtol主导可接受的误差。如果atol大于rtol * abs(y),则不能保证正确数字的数量。反之,为了达到期望的atol,设置rtol使得rtol * abs(y)始终小于atol可能是有益的。如果 y 的组成部分具有不同的比例,可能有益于通过传递形状为(n,)的 array_like 为atol的不同组件设置不同的atol值。rtol的默认值为 1e-3,atol的默认值为 1e-6。

vectorized布尔值,可选

fun是否可以以向量化方式调用。对于此求解器,建议设置为 False(默认)。

如果vectorized为 False,则fun始终使用形状为(n,)y调用,其中n = len(y0)

如果vectorized为 True,则fun可以使用形状为(n, k)y调用,其中k为整数。在这种情况下,fun必须表现出fun(t, y)[:, i] == fun(t, y[:, i])的行为(即返回数组的每一列是对应y列的状态的时间导数)。

设置vectorized=True允许‘Radau’和‘BDF’方法更快地近似雅可比矩阵的有限差分,但会导致此求解器执行速度较慢。

参考文献

[1]

P. Bogacki, L.F. Shampine, “A 3(2) Pair of Runge-Kutta Formulas”, Appl. Math. Lett. Vol. 2, No. 4. pp. 321-325, 1989.

属性:

nint

方程数量。

statusstring

求解器的当前状态:'running'(运行中)、'finished'(已完成)或 'failed'(失败)。

t_boundfloat

边界时间。

directionfloat

积分方向:+1 或 -1。

tfloat

当前时间。

yndarray

当前状态。

t_oldfloat

前一次时间。如果尚未进行步骤,则为 None。

step_sizefloat

最后一次成功步长的大小。如果尚未进行步骤,则为 None。

nfevint

系统右侧函数的评估次数。

njevint

雅可比矩阵的评估次数。对于此求解器始终为 0,因为它不使用雅可比矩阵。

nluint

LU 分解次数。对于此求解器始终为 0。

方法

dense_output() 计算在最后一次成功步骤上的局部插值。
step() 执行一步积分。

scipy.integrate.RK45

原文:docs.scipy.org/doc/scipy-1.12.0/reference/generated/scipy.integrate.RK45.html#scipy.integrate.RK45

class scipy.integrate.RK45(fun, t0, y0, t_bound, max_step=inf, rtol=0.001, atol=1e-06, vectorized=False, first_step=None, **extraneous)

显式五阶四阶 Runge-Kutta 方法。

这使用了 Dormand-Prince 一对公式 [1]。 错误控制假设精度为四阶方法的精度,但使用第五阶精确公式(进行局部外推)。 密集输出使用四次插值多项式 [2]

可应用于复域。

参数:

funcallable

系统的右手边。 调用签名为 fun(t, y)。 这里 t 是标量,ndarray y 有两种选择:它可以具有形状(n,);然后 fun 必须返回形状为(n,)的 array_like。 或者它可以具有形状(n,k);然后 fun 必须返回形状为(n,k)的 array_like,即每列对应于 y 中的单个列。 选择两种选项之间由 vectorized 参数决定(见下文)。

t0float

初始时间。

y0array_like,形状为(n,)

初始状态。

t_boundfloat

边界时间 - 积分不会超出它。 它还确定了积分的方向。

first_stepfloat 或 None,可选

初始步长。 默认为 None,这意味着算法应该选择。

max_stepfloat,可选

最大允许步长。 默认为 np.inf,即步长不受限制,完全由求解器确定。

rtol, atolfloat 和 array_like,可选

相对和绝对容差。 求解器使局部误差估计保持小于 atol + rtol * abs(y)。 这里 rtol 控制相对精度(正确数字的数量),而 atol 控制绝对精度(正确小数位的数量)。 要实现所需的 rtol,将 atol 设置为小于可以从 rtol * abs(y) 预期的最小值,以便 rtol 主导允许的误差。 如果 atol 大于 rtol * abs(y),则不能保证正确的数字。 相反,为了实现所需的 atol,设置 rtol,使得 rtol * abs(y) 始终小于 atol。 如果 y 的组件具有不同的尺度,则通过传递形状为(n,)的 array_like 来为 atol 的不同组件设置不同的 atol 值可能是有益的。 默认值为 rtol 为 1e-3 和 atol 为 1e-6。

vectorizedbool,可选

是否在矢量化方式中实现 fun。 默认为 False。

参考文献

[1]

J. R. Dormand, P. J. Prince,“A family of embedded Runge-Kutta formulae”,Journal of Computational and Applied Mathematics,Vol. 6,No. 1,pp. 19-26,1980。

[2]

L. W. Shampine,“Some Practical Runge-Kutta Formulas”,Mathematics of Computation,Vol. 46,No. 173,pp. 135-150,1986。

属性:

nint

方程的数量。

statusstring

求解器的当前状态:'running'、'finished'或'failed'。

t_boundfloat

边界时间。

directionfloat

积分方向:+1 或-1。

tfloat

当前时间。

yndarray

当前状态。

t_oldfloat

上一个时间。如果尚未进行任何步骤,则为 None。

step_sizefloat

最后一次成功步长的大小。如果尚未进行任何步骤,则为 None。

nfevint

系统右手边的评估次数。

njevint

雅可比矩阵的评估次数。对于这个求解器,始终为 0,因为它不使用雅可比矩阵。

nluint

LU 分解次数。对于这个求解器,始终为 0。

方法

dense_output() 计算在最后一次成功步骤上的局部插值。
step() 执行一次积分步骤。

scipy.integrate.DOP853

原文链接:docs.scipy.org/doc/scipy-1.12.0/reference/generated/scipy.integrate.DOP853.html#scipy.integrate.DOP853

class scipy.integrate.DOP853(fun, t0, y0, t_bound, max_step=inf, rtol=0.001, atol=1e-06, vectorized=False, first_step=None, **extraneous)

显式的 8 阶 Runge-Kutta 方法。

这是“DOP853”算法的 Python 实现,最初用 Fortran 编写[1][2]。请注意,这不是字面上的翻译,但算法核心和系数是相同的。

可以在复杂域中应用。

参数:

funcallable

系统的右侧。调用签名为fun(t, y)。这里,t是一个标量,而y是一个形状为(n,)的 ndarray 的两个选项之一:它可以返回形状为(n,)的 array_like,或者可以返回形状为(n, k)的 array_like,即每一列对应于y中的单个列。这两个选项的选择由vectorized参数决定(参见下文)。

t0float

初始时间。

y0array_like,形状为(n,)

初始状态。

t_boundfloat

边界时间 - 积分不会超出这个时间。它也决定了积分的方向。

first_stepfloat 或 None,可选

初始步长。默认为None,由算法选择。

max_stepfloat,可选

最大允许的步长。默认为 np.inf,即步长不受限制,完全由求解器确定。

rtol, atolfloat 和 array_like,可选

相对和绝对容差。求解器保持局部误差估计小于atol + rtol * abs(y)。这里,rtol控制相对精度(正确数字的数量),而atol控制绝对精度(正确小数位数)。为了达到期望的rtol,将atol设置为小于可以从rtol * abs(y)期望的最小值,以便rtol支配允许的误差。如果atol大于rtol * abs(y),则不能保证正确数字的数量。反之,为了达到期望的atol,设置rtol使得rtol * abs(y)始终小于atol可能是有益的。如果 y 的各个分量具有不同的比例,通过传递形状为(n,)的 array_like 的atol值为不同的分量设置不同的atol值。默认值为rtol为 1e-3 和atol为 1e-6。

vectorizedbool,可选

fun是否以向量化方式实现。默认值为 False。

参考文献

[1]

E. Hairer, S. P. Norsett G. Wanner,“求解普通微分方程 I:非刚性问题”,第 II 节。

[2]

DOP853 的原始 Fortran 代码页面

属性:

nint

方程的数量。

status字符串

求解器当前状态:‘running’,‘finished’或‘failed’。

t_boundfloat

边界时间。

directionfloat

积分方向:+1 或-1。

tfloat

当前时间。

yndarray

当前状态。

t_oldfloat

之前的时间。如果尚未进行步骤,则为无。

step_sizefloat

最后一个成功步骤的大小。如果尚未进行步骤,则为无。

nfevint

系统右侧的评估次数。

njevint

雅可比矩阵的评估次数。对于此求解器始终为 0,因为不使用雅可比矩阵。

nluint

LU 分解次数。对于此求解器始终为 0。

方法

dense_output() 计算最后一个成功步骤上的局部插值。
step() 执行一次积分步骤。

scipy.integrate.Radau

原文:docs.scipy.org/doc/scipy-1.12.0/reference/generated/scipy.integrate.Radau.html#scipy.integrate.Radau

class scipy.integrate.Radau(fun, t0, y0, t_bound, max_step=inf, rtol=0.001, atol=1e-06, jac=None, jac_sparsity=None, vectorized=False, first_step=None, **extraneous)

隐式龙格库塔方法,Radau IIA 5 阶族。

实现遵循 [1]。误差由三阶准确嵌入式公式控制。用于密集输出的满足配点条件的立方多项式。

参数:

funcallable

系统的右手边:状态 y 在时间 t 的时间导数。调用签名为 fun(t, y),其中 t 是标量,y 是形状为 len(y0) 的 ndarray。fun 必须返回与 y 相同形状的数组。有关更多信息,请参见 vectorized

t0float

初始时间。

y0array_like,形状为 (n,)

初始状态。

t_boundfloat

边界时间 - 集成不会超出此时间。它还决定了集成的方向。

first_stepfloat 或 None,可选

初始步长。默认为 None,表示算法应该选择。

max_stepfloat,可选

最大允许步长。默认为 np.inf,即步长不受限制,完全由求解器确定。

rtol, atolfloat 和 array_like,可选

相对和绝对容差。求解器保持局部误差估计小于 atol + rtol * abs(y)。这里 rtol 控制相对精度(正确数字的数量),而 atol 控制绝对精度(正确小数位的数量)。为了实现期望的 rtol,将 atol 设置为小于从 rtol * abs(y) 可预期的最小值,以使 rtol 主导允许的误差。如果 atol 大于 rtol * abs(y),则不能保证正确数字的数量。反之,为了实现期望的 atol,设置 rtol,使得 rtol * abs(y) 总是小于 atol。如果 y 的分量具有不同的比例,可能有利于通过传递形状为 (n,) 的 array_like 为 atol 的不同分量设置不同的 atol 值。默认值为 rtol 为 1e-3,atol 为 1e-6。

jac,可选

系统右手边的雅可比矩阵相对于 y,该方法所需。雅可比矩阵的形状为 (n, n),其元素 (i, j) 等于 d f_i / d y_j。有三种定义雅可比矩阵的方法:

  • 如果是 array_like 或 sparse_matrix,则假定雅可比矩阵是恒定的。
  • 如果为 callable,则假定雅可比矩阵依赖于 t 和 y;将按需调用为 jac(t, y)。对于 'Radau' 和 'BDF' 方法,返回值可能是稀疏矩阵。
  • 如果为 None(默认),则雅可比矩阵将通过有限差分近似。

通常建议提供雅可比矩阵,而不是依赖有限差分近似。

jac_sparsity,可选

为有限差分近似的雅可比矩阵定义稀疏结构。其形状必须为(n, n)。如果jac不是None,则忽略此参数。如果雅可比矩阵在每行中只有少量非零元素,提供稀疏结构将极大地加快计算速度[2]。零条目表示雅可比矩阵中相应元素始终为零。如果为 None(默认),则假定雅可比矩阵是密集的。

vectorizedbool,可选

fun是否可以以向量化方式调用的能力。默认为 False。

如果vectorized为 False,则fun将始终使用形状为(n,)y调用,其中n = len(y0)

如果vectorized为 True,则可以用形状为(n, k)y调用fun,其中k是整数。在这种情况下,fun必须表现出fun(t, y)[:, i] == fun(t, y[:, i])(即返回数组的每一列都是与y的一列对应的状态的时间导数)。

设置vectorized=True允许通过此方法更快地进行雅可比矩阵的有限差分近似,但在某些情况下可能会导致整体执行速度较慢(例如,len(y0)较小)。

参考文献

[1]

E. Hairer, G. Wanner,“解常微分方程 II:刚性和微分代数问题”,第 IV.8 节。

[2]

A. Curtis, M. J. D. Powell, and J. Reid,“关于稀疏雅可比矩阵估计的研究”,应用数学研究院杂志,13,pp. 117-120,1974。

属性:

nint

方程数量。

statusstring

求解器的当前状态:‘running’、‘finished’或‘failed’。

t_boundfloat

边界时间。

directionfloat

积分方向:+1 或-1。

tfloat

当前时间。

yndarray

当前状态。

t_oldfloat

上一个时间。如果尚未执行任何步骤,则为 None。

step_sizefloat

最后一次成功步长的大小。如果尚未执行任何步骤,则为 None。

nfevint

右手边函数评估次数。

njevint

雅可比矩阵的评估次数。

nluint

LU 分解次数。

方法

dense_output() 计算上一次成功步骤的局部插值。
step() 执行一次积分步骤。

scipy.integrate.BDF

原文:docs.scipy.org/doc/scipy-1.12.0/reference/generated/scipy.integrate.BDF.html#scipy.integrate.BDF

class scipy.integrate.BDF(fun, t0, y0, t_bound, max_step=inf, rtol=0.001, atol=1e-06, jac=None, jac_sparsity=None, vectorized=False, first_step=None, **extraneous)

基于向后差分公式的隐式方法。

这是一个变阶方法,阶数自动从 1 变化到 5。BDF 算法的一般框架描述在[1]中。该类实现了一种准恒定步长,如[2]所述。常步长 BDF 的误差估计策略在[3]中推导。还实现了使用修改的公式(NDF)增强精度[2]

可应用于复杂域。

参数:

funcallable

系统的右手边:状态y在时间t的时间导数。调用签名是fun(t, y),其中t是标量,y是形状为len(y0)的 ndarray。fun必须返回与y相同形状的数组。详见向量化获取更多信息。

t0float

初始时间。

y0array_like,形状为(n,)

初始状态。

t_boundfloat

边界时间 - 积分不会超出此时间。它还决定了积分的方向。

first_stepfloat 或 None,可选

初始步长。默认为None,表示算法应选择。

max_stepfloat, 可选

最大允许步长。默认为 np.inf,即步长不受限制,完全由求解器决定。

rtol, atolfloat 和 array_like,可选

相对和绝对容差。求解器保持本地误差估计小于atol + rtol * abs(y)。这里rtol控制相对精度(正确位数),而atol控制绝对精度(正确小数位数)。为了达到期望的rtol,将atol设为比从rtol * abs(y)预期的最小值小,使得rtol主导可接受的误差。如果atol大于rtol * abs(y),则不能保证正确位数。相反,为了达到期望的atol,设置rtol,使得rtol * abs(y)始终小于atol。如果 y 的分量具有不同的比例,通过传递形状为(n,)的 array_like 给atol,为不同的分量设置不同的atol值可能是有益的。默认值为 1e-3(rtol)和 1e-6(atol)。

jac,可选

系统右侧的雅可比矩阵与 y 的关系,该方法所需。雅可比矩阵形状为(n, n),其元素(i, j)等于d f_i / d y_j。有三种定义雅可比矩阵的方法:

  • 如果是 array_like 或 sparse_matrix,则假定雅可比矩阵是常数。
  • 如果可调用,则假定雅可比矩阵依赖于 t 和 y;将根据需要调用为jac(t, y)。对于‘Radau’和‘BDF’方法,返回值可能是稀疏矩阵。
  • 如果为 None(默认),雅可比矩阵将通过有限差分逼近来近似。

通常建议提供雅可比矩阵,而不是依赖有限差分逼近。

jac_sparsity,可选

为有限差分逼近的雅可比矩阵定义稀疏结构。其形状必须为(n, n)。如果jac不是None,则此参数将被忽略。如果雅可比矩阵每行只有少数非零元素,则提供稀疏结构将极大地加速计算 [4]。零项表示雅可比矩阵中的对应元素始终为零。如果为 None(默认),则假定雅可比矩阵为密集型。

vectorizedbool,可选

fun是否可以以矢量化方式调用的标志。默认为 False。

如果vectorized为 False,fun将始终使用形状为(n,)y调用,其中n = len(y0)

如果vectorized为 True,则fun可以使用形状为(n, k)y调用,其中k是整数。在这种情况下,fun必须表现出fun(t, y)[:, i] == fun(t, y[:, i])(即返回数组的每一列都是与y的每一列对应的状态的时间导数)。

设置vectorized=True允许通过此方法更快地进行雅可比矩阵的有限差分逼近,但在某些情况下(例如len(y0)较小)可能导致总体执行速度较慢。

参考文献

[1]

G. D. Byrne, A. C. Hindmarsh,“用于数值解普通微分方程的多算法”,ACM Transactions on Mathematical Software,Vol. 1,No. 1,pp. 71-96,1975 年 3 月。

[2] (1,2)

L. F. Shampine, M. W. Reichelt,“MATLAB ODE SUITE”,SIAM J. SCI. COMPUTE.,Vol. 18,No. 1,pp. 1-22,1997 年 1 月。

[3]

E. Hairer, G. Wanner,“求解普通微分方程 I:非刚性问题”,第 III.2 节。

[4]

A. Curtis, M. J. D. Powell, 和 J. Reid,“关于稀疏雅可比矩阵估计的问题”,数学应用研究所学报,13,pp. 117-120,1974。

属性:

nint

方程数量。

statusstring

求解器的当前状态:‘running’、‘finished’或‘failed’。

t_boundfloat

边界时间。

directionfloat

积分方向:+1 或-1。

tfloat

当前时间。

yndarray

当前状态。

t_oldfloat

上一个时间。如果尚未进行步骤,则为 None。

step_sizefloat

最后一个成功步长的大小。如果尚未进行步骤,则为 None。

nfevint

右侧函数评估次数。

njevint

雅可比矩阵的评估次数。

nluint

LU 分解次数。

方法

dense_output() 计算最后一个成功步骤上的本地插值器。
step() 执行一步积分。

scipy.integrate.LSODA

原文链接:docs.scipy.org/doc/scipy-1.12.0/reference/generated/scipy.integrate.LSODA.html#scipy.integrate.LSODA

class scipy.integrate.LSODA(fun, t0, y0, t_bound, first_step=None, min_step=0.0, max_step=inf, rtol=0.001, atol=1e-06, jac=None, lband=None, uband=None, vectorized=False, **extraneous)

带有自动刚性检测和切换的 Adams/BDF 方法。

这是一个对来自 ODEPACK 的 Fortran 求解器的包装器 [1]。它在非刚性 Adams 方法和刚性 BDF 方法之间自动切换。该方法最初在 [2] 中详细描述。

参数:

funcallable

系统的右手边:时间 t 处状态 y 的时间导数。调用签名为 fun(t, y),其中 t 是标量,y 是形状为 len(y0) 的 ndarray。fun 必须返回与 y 相同形状的数组。更多信息请参见向量化

t0float

初始时间。

y0array_like,形状为 (n,)

初始状态。

t_boundfloat

边界时间 - 积分不会超出此时间。它还决定了积分的方向。

first_stepfloat 或 None,可选

初始步长。默认为 None,表示算法应选择。

min_stepfloat,可选

允许的最小步长。默认为 0.0,即步长不受限制,完全由求解器确定。

max_stepfloat,可选

允许的最大步长。默认为 np.inf,即步长不受限制,完全由求解器确定。

rtol, atolfloat 和 array_like,可选

相对和绝对容差。求解器保持局部误差估计小于 atol + rtol * abs(y)。这里 rtol 控制相对精度(正确数字的数量),而 atol 控制绝对精度(正确小数位数)。为了实现期望的 rtol,设置 atol 小于从 rtol * abs(y) 可预期的最小值,以便 rtol 主导可允许的误差。如果 atol 大于 rtol * abs(y),则不能保证正确数字的数量。相反,为了实现期望的 atol,设置 rtol,使得 rtol * abs(y) 总是小于 atol。如果 y 的各分量具有不同的尺度,则通过传递形状为 (n,) 的 array_like 来为 atol 的不同分量设置不同的值可能是有益的。默认值为 rtol 的 1e-3 和 atol 的 1e-6。

jacNone 或 callable,可选

系统右手边关于 y 的雅可比矩阵。雅可比矩阵形状为 (n, n),其元素 (i, j) 等于 d f_i / d y_j。函数将作为 jac(t, y) 调用。如果为 None(默认),雅可比将通过有限差分近似。通常建议提供雅可比矩阵,而不是依赖于有限差分近似。

lband, ubandint 或 None

定义雅可比矩阵带宽的参数,即,jac[i, j] != 0 仅当 i - lband <= j <= i + uband。设置这些参数要求您的 jac 程序以压缩格式返回雅可比矩阵:返回的数组必须具有 n 列和 uband + lband + 1 行,其中雅可比矩阵的对角线被写入。具体而言,jac_packed[uband + i - j , j] = jac[i, j]。同样的格式也用于 scipy.linalg.solve_banded(请参考示例)。这些参数也可以与 jac=None 一起使用,以减少通过有限差分估计的雅可比元素数量。

vectorized bool,可选

fun 是否可以以矢量化方式调用。建议此求解器默认为 False。

如果 vectorized 为 False,则 fun 将始终使用形状为 (n,)y 调用,其中 n = len(y0)

如果 vectorized 为 True,则 fun 可能以形状为 (n, k)y 调用,其中 k 是整数。在这种情况下,fun 必须使得 fun(t, y)[:, i] == fun(t, y[:, i])(即返回数组的每一列都是与 y 的相应列对应的状态的时间导数)。

设置 vectorized=True 允许 ‘Radau’ 和 ‘BDF’ 方法更快地通过有限差分逼近雅可比矩阵,但会导致此求解器执行速度较慢。

参考文献

[1]

A. C. Hindmarsh, “ODEPACK, A Systematized Collection of ODE Solvers,” IMACS Transactions on Scientific Computation, Vol 1., pp. 55-64, 1983.

[2]

L. Petzold, “Automatic selection of methods for solving stiff and nonstiff systems of ordinary differential equations”, SIAM Journal on Scientific and Statistical Computing, Vol. 4, No. 1, pp. 136-148, 1983.

属性:

n int

方程的数量。

status string

求解器的当前状态:‘running’(运行中)、‘finished’(已完成)或‘failed’(失败)。

t_bound float

边界时间。

direction float

积分方向:+1 或 -1。

t float

当前时间。

y ndarray

当前状态。

t_old float

前一个时间。如果还没有进行步骤,则为无。

nfev int

右侧求值的次数。

njev int

雅可比矩阵的求值次数。

方法

dense_output() 计算上一次成功步骤的局部插值。
step() 执行一步积分。

scipy.integrate.OdeSolver

原文链接:docs.scipy.org/doc/scipy-1.12.0/reference/generated/scipy.integrate.OdeSolver.html#scipy.integrate.OdeSolver

class scipy.integrate.OdeSolver(fun, t0, y0, t_bound, vectorized, support_complex=False)

ODE 求解器的基类。

要实现新的求解器,需要遵循以下准则:

  1. 构造函数必须接受在基类中呈现的参数(下面列出),以及与求解器特定的任何其他参数。
  2. 构造函数必须接受任意多余参数**extraneous,但通过common.warn_extraneous函数警告这些参数是不相关的。不要将这些参数传递给基类。
  3. 求解器必须实现一个私有方法_step_impl(self),将求解器推进一步。必须返回元组(success, message),其中success是一个布尔值,指示步骤是否成功,message是包含失败描述的字符串(如果步骤失败)或 None。
  4. 求解器必须实现一个私有方法_dense_output_impl(self),返回一个DenseOutput对象,覆盖最后一个成功步骤。
  5. 求解器必须具有以下属性列表中列出的属性。注意,t_oldstep_size会自动更新。
  6. 使用fun(self, t, y)方法进行系统右手边评估,这样函数评估数(nfev)会自动跟踪。
  7. 为方便起见,基类提供了fun_single(self, t, y)fun_vectorized(self, t, y),分别用于非向量化和向量化方式评估右手边(不管构造函数中的fun如何实现)。这些调用不会增加nfev
  8. 如果求解器使用雅可比矩阵和 LU 分解,它应该追踪雅可比矩阵评估数(njev)和 LU 分解数(nlu)。
  9. 根据惯例,用于计算雅可比矩阵有限差分近似的函数评估不应计入nfev,因此在计算雅可比矩阵有限差分近似时,请使用fun_single(self, t, y)fun_vectorized(self, t, y)

参数:

funcallable

系统右手边:时间t处状态y的时间导数。调用签名为fun(t, y),其中t是标量,y是具有len(y) = len(y0)的 ndarray。fun必须返回与y相同形状的数组。有关更多信息,请参见vectorized

t0float

初始时间。

y0array_like,形状为(n,)

初始状态。

t_boundfloat

边界时间 —— 积分不会超出它。它还确定积分的方向。

vectorizedbool

fun是否可以以向量化方式调用。默认为 False。

如果vectorized为 False,则fun始终以形状为(n,)y调用,其中n = len(y0)

如果vectorized为 True,则可以使用形状为(n, k)y调用fun,其中k是整数。在这种情况下,fun必须表现出fun(t, y)[:, i] == fun(t, y[:, i])(即返回数组的每一列是与y的一列对应的状态的时间导数)。

设置vectorized=True允许方法‘Radau’和‘BDF’通过更快的有限差分逼近雅可比矩阵,但会导致其他方法执行较慢。在某些情况下(例如y0很小),它也可能导致‘Radau’和‘BDF’的整体执行较慢。

support_complex 布尔值,可选

是否应支持复数域中的积分。通常由派生的求解器类能力决定。默认为 False。

属性:

n 整数

方程的数量。

status 字符串

求解器的当前状态:‘运行中’,‘已完成’或‘失败’。

t_bound 浮点数

边界时间。

direction 浮点数

积分方向:+1 或 -1。

t 浮点数

当前时间。

y 数组

当前状态。

t_old 浮点数

先前时间。如果尚未执行步骤,则为 None。

step_size 浮点数

上一个成功步骤的大小。如果尚未执行步骤,则为 None。

nfev 整数

系统右手边评估的数量。

njev 整数

雅可比矩阵评估的数量。

nlu 整数

LU 分解的数量。

方法

dense_output() 计算上一次成功步骤的局部插值。
step() 执行一步积分。

scipy.integrate.DenseOutput

原文链接:docs.scipy.org/doc/scipy-1.12.0/reference/generated/scipy.integrate.DenseOutput.html#scipy.integrate.DenseOutput

class scipy.integrate.DenseOutput(t_old, t)

ODE 求解器生成的局部插值器的基类。

它在t_mint_max之间进行插值(见下面的属性)。超出此区间的评估不被禁止,但精度不能保证。

属性:

t_min, t_maxfloat

插值的时间范围。

方法

__call__(t) 评估插值函数。

scipy.integrate.OdeSolution

原文:docs.scipy.org/doc/scipy-1.12.0/reference/generated/scipy.integrate.OdeSolution.html#scipy.integrate.OdeSolution

class scipy.integrate.OdeSolution(ts, interpolants, alt_segment=False)

连续 ODE 解决方案。

它组织为一组DenseOutput对象,代表局部插值器。 它提供了一个算法来为每个给定点选择合适的插值器。

插值器覆盖从t_mint_max的范围(见下面的属性)。 虽然不禁止在此间隔之外进行评估,但不能保证准确性。

在断点(ts中的一个值)处评估时,将选择具有较低索引的段。

参数:

tsarray_like,形状为(n_segments + 1,)

定义局部插值器的时间点。 必须严格递增或递减(允许两点的零段)。

interpolantsDenseOutput 对象列表,具有 n_segments 个元素

局部插值器。 假定第 i 个插值器在ts[i]ts[i + 1]之间定义。

alt_segment布尔值

请求备选插值器段选择方案。 在每个求解器积分点上,两个插值器段可用。 默认(False)和备选(True)行为分别选择所请求时间对应的段与t_old。 此功能仅适用于测试插值器的准确性:不同的积分器使用不同的构造策略。

属性:

t_min, t_max浮点数

插值的时间范围。

方法

__call__(t) 评估解决方案。

scipy.integrate.odeint

原文:docs.scipy.org/doc/scipy-1.12.0/reference/generated/scipy.integrate.odeint.html#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, …) 或 callable(t, y, …)

在 t 处计算 y 的导数。如果签名是 callable(t, y, ...),则参数 tfirst 必须设置为 True

y0数组

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

t数组

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

args元组,可选

传递给函数的额外参数。

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

func 的梯度(雅可比矩阵)。如果签名是 callable(t, y, ...),则参数 tfirst 必须设置为 True

col_derivbool,可选

如果 Dfun 定义沿列的导数(更快),否则 Dfun 应定义沿行的导数。

full_outputbool,可选

如果返回一个字典作为第二个输出的可选输出,则为真

printmessgbool,可选

是否打印收敛消息

tfirstbool,可选

如果为真,则 func(和 Dfun(如果给定))的前两个参数必须为 t, y,而不是默认的 y, t

新版本 1.1.0 中的新增功能。

返回:

y数组,形状为 (len(t), len(y0))

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

infodictdict,仅在 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应返回一个矩阵,其行包含非零带(从最低对角线开始)。因此,来自Dfun的返回矩阵jac应具有形状(ml + mu + 1, len(y0)),当ml >=0mu >=0时。 jac中的数据必须存储,以便jac[i - j + mu, j]保存第i个方程相对于第j个状态变量的导数。如果col_deriv为 True,则必须返回此jac的转置。

rtol, atolfloat, optional

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

tcritndarray, optional

关键点(例如奇点)的向量,需要对积分进行注意。

h0float, (0: solver-determined), optional

尝试在第一步上尝试的步长。

hmaxfloat, (0: solver-determined), optional

允许的最大绝对步长大小。

hminfloat, (0: solver-determined), optional

允许的最小绝对步长大小。

ixprbool, optional

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

mxstepint, (0: solver-determined), optional

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

mxhnilint, (0: solver-determined), optional

允许打印的最大消息数量。

mxordnint, (0: solver-determined), optional

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

mxordsint, (0: solver-determined), optional

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

另见

solve_ivp

解决 ODE 系统的初始值问题

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

scipy.integrate.ode

原文链接:docs.scipy.org/doc/scipy-1.12.0/reference/generated/scipy.integrate.ode.html#scipy.integrate.ode

class scipy.integrate.ode(f, jac=None)

一个通用的数值积分器接口类。

解一个方程系统 (y'(t) = f(t,y)) ,可选参数 jac = df/dy

f(t, y, ...) 的前两个参数与 scipy.integrate.odeint 中系统定义函数中的参数顺序相反。

参数:

fcallable f(t, y, *f_args)

微分方程的右侧。t 是一个标量,y.shape == (n,)。通过调用 set_f_params(*args) 来设置 f_argsf 应返回标量、数组或列表(而不是元组)。

jaccallable jac(t, y, *jac_args),可选

右侧的雅可比矩阵,jac[i,j] = d f[i] / d y[j]。通过调用 set_jac_params(*args) 来设置 jac_args

另请参见

odeint

一个基于 ODEPACK 中 lsoda 的简单接口的积分器。

quad

用于求曲线下面积的工具

注意

可用的积分器如下所示。可以使用 set_integrator 方法选择它们。

“vode”

实数变系数普通微分方程求解器,具有固定前导系数实现。它提供了隐式的亚当斯方法(用于非刚性问题)和基于向后差分公式(BDF)的方法(用于刚性问题)。

来源:www.netlib.org/ode/vode.f

警告

该积分器不可重入。你不能同时使用两个使用“vode”积分器的 ode 实例。

该积分器在 ode 类的 set_integrator 方法中接受以下参数:

  • atol:float 或 sequence 解的绝对容差
  • rtol:float 或 sequence 解的相对容差
  • lband:None 或 int
  • uband:None 或 int 雅可比带宽,对于 i-lband <= j <= i+uband,jac[i,j] != 0。设置这些需要你的 jac 程序以打包格式返回雅可比矩阵,jac_packed[i-j+uband, j] = jac[i,j]。矩阵的维度必须是 (lband+uband+1, len(y))。
  • method: ‘adams’ 或 ‘bdf’ 选择使用的求解器,Adams(非刚性)或 BDF(刚性)
  • with_jacobian : bool 此选项仅在用户未提供雅可比函数且未指示(通过设置任何带状)雅可比矩阵为带状时考虑。在这种情况下,with_jacobian指定了 ODE 求解器校正步骤的迭代方法,可以是使用内部生成的完整雅可比矩阵的弦迭代,或者是不使用雅可比矩阵的功能迭代。
  • nsteps : int 一次调用求解器期间允许的最大(内部定义的)步数。
  • first_step : float
  • min_step : float
  • max_step : float 积分器使用的步长限制。
  • order : int 积分器使用的最大阶数,对于 Adams 阶数 <= 12,对于 BDF 阶数 <= 5。

“zvode”

复值变系数常微分方程求解器,带有固定前导系数实现。它提供隐式的 Adams 方法(用于非刚性问题)和基于后向差分公式(BDF)的方法(用于刚性问题)。

Source: www.netlib.org/ode/zvode.f

警告

此积分器不是可重入的。您不能同时使用两个 ode 实例来使用“zvode”积分器。

此积分器在 set_integrator 中接受与“vode”求解器相同的参数。

注意

当在刚性系统中使用 ZVODE 时,应仅用于函数 f 是解析的情况,即每个 f(i)是每个 y(j)的解析函数。解析性意味着偏导数 df(i)/dy(j)是唯一的复数,并且这一事实对 ZVODE 解决刚性情况下出现的密集或带状线性系统至关重要。对于一个复杂的刚性 ODE 系统,其中 f 不是解析的情况,ZVODE 可能会出现收敛失败,对于这个问题,应该使用等效的实系统中的 DVODE。

“lsoda”

实值变系数常微分方程求解器,带有固定前导系数实现。它提供在隐式 Adams 方法(用于非刚性问题)和基于后向差分公式(BDF)的方法(用于刚性问题)之间的自动方法切换。

Source: www.netlib.org/odepack

警告

此积分器不是可重入的。您不能同时使用两个 ode 实例来使用“lsoda”积分器。

此积分器在 set_integrator 方法中的 ode 类中接受以下参数:

  • atol : float 或序列 解的绝对容差
  • rtol : float 或序列 解的相对容差
  • lband : None 或 int
  • uband : None 或 int 雅可比矩阵的带宽,jac[i,j] != 0 对于 i-lband <= j <= i+uband。设置这些需要您的 jac 例程以紧凑格式返回雅可比矩阵,jac_packed[i-j+uband, j] = jac[i,j]。
  • with_jacobian : bool 未使用。
  • nsteps : int 在一次调用解算器期间允许的最大(内部定义的)步数。
  • first_step : float
  • min_step : float
  • max_step : float 集成器使用的步长限制。
  • max_order_ns : int 在非刚性情况下使用的最大阶数(默认为 12)。
  • max_order_s : int 在刚性情况下使用的最大阶数(默认为 5)。
  • max_hnil : int 报告步长过小的消息数的最大数目(t + h = t)(默认为 0)
  • ixpr : int 是否在方法切换时生成额外的打印输出(默认为 False)。

“dopri5”

这是一种显式 Runge-Kutta 方法,阶数为(4)5,由 Dormand 和 Prince 提出(具有步长控制和密集输出)。

作者:

E. Hairer 和 G. Wanner 瑞士日内瓦大学,数学系 CH-1211 Geneve 24,瑞士 电子邮件:ernst.hairer@math.unige.ch,gerhard.wanner@math.unige.ch

本代码在[HNW93]中有描述。

此集成器在 ode 类的 set_integrator()方法中接受以下参数:

  • atol : float 或序列的解的绝对容差
  • rtol : float 或序列的解的相对容差
  • nsteps : int 在一次调用解算器期间允许的最大(内部定义的)步数。
  • first_step : float
  • max_step : float
  • safety : float 对新步长选择的安全因子(默认为 0.9)
  • ifactor : float
  • dfactor : float 在一个步骤中增加/减少步长的最大因子。
  • beta : float 控制稳定步长的 Beta 参数。
  • verbosity : int 用于打印消息的开关(小于 0 表示不打印消息)。

“dop853”

这是一种由 Dormand 和 Prince 提出的显式 Runge-Kutta 方法,阶数为 8(5,3)(具有步长控制和密集输出)。

选项和引用与“dopri5”相同。

参考文献

[HNW93]

E. Hairer, S.P. Norsett 和 G. Wanner,《求解常微分方程》第二版。Springer 计算数学系列,Springer-Verlag(1993 年)

示例

一个集成问题及其相应的雅可比矩阵:

>>> from scipy.integrate import ode
>>>
>>> y0, t0 = [1.0j, 2.0], 0
>>>
>>> def f(t, y, arg1):
...     return [1j*arg1*y[0] + y[1], -arg1*y[1]**2]
>>> def jac(t, y, arg1):
...     return [[1j*arg1, 1], [0, -arg1*2*y[1]]] 

集成:

>>> r = ode(f, jac).set_integrator('zvode', method='bdf')
>>> r.set_initial_value(y0, t0).set_f_params(2.0).set_jac_params(2.0)
>>> t1 = 10
>>> dt = 1
>>> while r.successful() and r.t < t1:
...     print(r.t+dt, r.integrate(r.t+dt))
1 [-0.71038232+0.23749653j  0.40000271+0.j        ]
2.0 [0.19098503-0.52359246j 0.22222356+0.j        ]
3.0 [0.47153208+0.52701229j 0.15384681+0.j        ]
4.0 [-0.61905937+0.30726255j  0.11764744+0.j        ]
5.0 [0.02340997-0.61418799j 0.09523835+0.j        ]
6.0 [0.58643071+0.339819j 0.08000018+0.j      ]
7.0 [-0.52070105+0.44525141j  0.06896565+0.j        ]
8.0 [-0.15986733-0.61234476j  0.06060616+0.j        ]
9.0 [0.64850462+0.15048982j 0.05405414+0.j        ]
10.0 [-0.38404699+0.56382299j  0.04878055+0.j        ] 

属性:

tfloat

当前时间。

yndarray

当前变量值。

方法

get_return_code() 提取集成的返回代码,以便在集成失败时进行更好的控制。
integrate(t[, step, relax]) 找到 y=y(t),将 y 设置为初始条件,并返回 y。
set_f_params(*args) 为用户提供的函数 f 设置额外的参数。
set_initial_value(y[, t]) 设置初始条件 y(t) = y。
set_integrator(name, **integrator_params) 根据名称设置积分器。
set_jac_params(*args) 为用户提供的函数 jac 设置额外参数。
set_solout(solout) 设置在每次成功积分步骤时调用的可调用函数。
successful() 检查积分是否成功。
posted @ 2024-06-27 17:11  绝不原创的飞龙  阅读(88)  评论(0)    收藏  举报