SciPy 1.12 中文文档（十）

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

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{None, array_like, sparse_matrix, callable}，可选

系统右侧的雅可比矩阵与 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{None, array_like, 稀疏矩阵}，可选

为有限差分逼近的雅可比矩阵定义稀疏结构。其形状必须为(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. 构造函数必须接受在基类中呈现的参数（下面列出），以及与求解器特定的任何其他参数。
3. 构造函数必须接受任意多余参数**extraneous，但通过common.warn_extraneous函数警告这些参数是不相关的。不要将这些参数传递给基类。
5. 求解器必须实现一个私有方法*_step_impl(self)*，将求解器推进一步。必须返回元组(success, message)，其中success是一个布尔值，指示步骤是否成功，message是包含失败描述的字符串（如果步骤失败）或 None。
7. 求解器必须实现一个私有方法*_dense_output_impl(self)*，返回一个DenseOutput对象，覆盖最后一个成功步骤。
9. 求解器必须具有以下属性列表中列出的属性。注意，t_old和step_size会自动更新。
11. 使用*fun(self, t, y)*方法进行系统右手边评估，这样函数评估数(nfev)会自动跟踪。
13. 为方便起见，基类提供了fun_single(self, t, y)和fun_vectorized(self, t, y)，分别用于非向量化和向量化方式评估右手边（不管构造函数中的fun如何实现）。这些调用不会增加nfev。
15. 如果求解器使用雅可比矩阵和 LU 分解，它应该追踪雅可比矩阵评估数(njev)和 LU 分解数(nlu)。
17. 根据惯例，用于计算雅可比矩阵有限差分近似的函数评估不应计入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_min和t_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_min到t_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 >=0或mu >=0时。 jac中的数据必须存储，以便jac[i - j + mu, j]保存第i个方程相对于第j个状态变量的导数。如果col_deriv为 True，则必须返回此jac的转置。

rtol, atolfloat, optional

输入参数rtol和atol确定求解器执行的误差控制。求解器将根据形式为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 

其中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() 

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_args。f 应返回标量、数组或列表（而不是元组）。

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()	检查积分是否成功。

scipy.integrate.complex_ode

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

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

用于复杂系统的 ode 的包装器。

这个函数类似于ode，但在使用积分器之前，将复值方程系统重新映射为实值方程系统。

参数：

f可调用函数 f(t, y, *f_args)

方程的右手边。t 是标量，y.shape == (n,)。通过调用 set_f_params(*args) 设置 f_args。

jac可调用函数 jac(t, y, *jac_args)

方程的雅可比矩阵，jac[i,j] = d f[i] / d y[j]。通过调用 set_f_params(*args) 设置 jac_args。

示例

有关用法示例，请参见ode。

属性：

t浮点数

当前时间。

y数组

当前变量值。

方法

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()	检查积分是否成功。

scipy.integrate.ODEintWarning

原文：docs.scipy.org/doc/scipy-1.12.0/reference/generated/scipy.integrate.ODEintWarning.html#scipy.integrate.ODEintWarning

exception scipy.integrate.ODEintWarning

在执行odeint时引发警告。

with_traceback()

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

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 解决一组普通微分方程。

解决第一阶 ODE-s 的初始值问题，对于刚性或非刚性系统：

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_deriv布尔值，可选

如果Dfun沿列定义导数（更快），否则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, optional

如果这两者都不是 None 或非负，则假定雅可比矩阵为带状矩阵。这些数字给出此带状矩阵中的下限和上限非零对角线数。对于带状情况，Dfun应返回一个矩阵，其行包含非零带（从最低对角线开始）。因此，来自Dfun的返回矩阵jac在ml >=0或mu >=0时应具有形状(ml + mu + 1, len(y0))。jac中的数据必须存储为jac[i - j + mu, j]，其中i方程对j状态变量的导数。如果col_deriv为 True，则必须返回此jac的转置。

rtol, atolfloat, optional

输入参数rtol和atol确定求解器执行的误差控制。求解器将根据形式为max-norm of (e / ewt) <= 1的不等式控制 y 的估计局部误差向量 e，其中 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

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

[quad`

用于找到曲线下面积

示例

受重力和摩擦影响的钟摆角度 theta 的二阶微分方程可以写成：

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

其中 b 和 c 是正常数，prime（’）表示导数。要用 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() 

scipy.integrate.solve_bvp

原文：docs.scipy.org/doc/scipy-1.12.0/reference/generated/scipy.integrate.solve_bvp.html#scipy.integrate.solve_bvp

scipy.integrate.solve_bvp(fun, bc, x, y, p=None, S=None, fun_jac=None, bc_jac=None, tol=0.001, max_nodes=1000, verbose=0, bc_tol=None)

解决 ODE 系统的边界值问题。

此函数数值解一个带有两点边界条件的一阶 ODE 系统：

dy / dx = f(x, y, p) + S * y / (x - a), a <= x <= b
bc(y(a), y(b), p) = 0 

这里 x 是一个 1-D 独立变量，y(x) 是一个 N-D 向量值函数，p 是一个 k-D 未知参数向量，它与 y(x) 一起被找到。为了确定问题，必须有 n + k 个边界条件，即 bc 必须是一个 (n + k)-D 函数。

系统右手边的最后一个奇异项是可选的。它由一个 n×n 矩阵 S 定义，使得解必须满足 S y(a) = 0。此条件将在迭代过程中强制执行，因此不得与边界条件相矛盾。详见 [2]，解释在数值求解 BVPs 时如何处理此项。

也可以解决复杂域中的问题。在这种情况下，y 和 p 被视为复数，f 和 bc 被假定为复值函数，但 x 保持实数。注意 f 和 bc 必须是复可微的（满足柯西-黎曼方程 [4]），否则应将问题分别重写为实部和虚部。要在复杂域中解决问题，请传递一个带有复数数据类型的初始猜测值 y。

参数：

funcallable

系统的右手边。调用签名为 fun(x, y) 或者如果存在参数则为 fun(x, y, p)。所有参数都是 ndarray：x 的形状为 (m,)，y 的形状为 (n, m)，意味着 y[:, i] 对应于 x[i]，p 的形状为 (k,)。返回值必须是形状为 (n, m) 的数组，并且与 y 的布局相同。

bccallable

评估边界条件残差的函数。调用签名为 bc(ya, yb) 或者如果存在参数则为 bc(ya, yb, p)。所有参数都是 ndarray：ya 和 yb 的形状为 (n,)，p 的形状为 (k,)。返回值必须是形状为 (n + k,) 的数组。

xarray_like，形状为 (m,)

初始网格。必须是一系列严格增加的实数，满足 x[0]=a 和 x[-1]=b。

yarray_like，形状为 (n, m)

函数在网格节点处的初始猜测值，第 i 列对应于 x[i]。对于复数域中的问题，即使初始猜测是纯实数，也要传递带有复数数据类型的 y。

p形状为 (k,) 的数组或 None，可选

未知参数的初始猜测值。如果为 None（默认），则假定问题不依赖于任何参数。

S形状为 (n, n) 的数组或 None

定义奇异项的矩阵。如果为 None（默认），则在没有奇异项的情况下解决问题。

fun_jaccallable 或 None，可选

计算 f 对 y 和 p 的导数的函数。其调用签名为 fun_jac(x, y) 或者如果存在参数则为 fun_jac(x, y, p)。返回必须按以下顺序包含 1 或 2 个元素：

• df_dy：形状为 (n, n, m) 的 array_like，其中元素 (i, j, q) 等于 d f_i(x_q, y_q, p) / d (y_q)_j。
• df_dp：形状为 (n, k, m) 的 array_like，其中元素 (i, j, q) 等于 d f_i(x_q, y_q, p) / d p_j。

此处 q 表示 x 和 y 定义的节点数，而 i 和 j 表示向量分量数。如果问题在没有未知参数的情况下解决，则不应返回 df_dp。

如果 fun_jac 为 None（默认情况下），则通过向前有限差分法估计导数。

bc_jac可调用对象或 None，可选

计算 bc 对 ya、yb 和 p 的导数的函数。其调用签名为 bc_jac(ya, yb) 或者如果存在参数则为 bc_jac(ya, yb, p)。返回必须按以下顺序包含 2 或 3 个元素：

• dbc_dya：形状为 (n, n) 的 array_like，其中元素 (i, j) 等于 d bc_i(ya, yb, p) / d ya_j。
• dbc_dyb：形状为 (n, n) 的 array_like，其中元素 (i, j) 等于 d bc_i(ya, yb, p) / d yb_j。
• dbc_dp：形状为 (n, k) 的 array_like，其中元素 (i, j) 等于 d bc_i(ya, yb, p) / d p_j。

如果问题在没有未知参数的情况下解决，则不应返回 dbc_dp。

如果 bc_jac 为 None（默认情况下），则通过向前有限差分法估计导数。

tolfloat，可选

求解的所需解的容差。如果我们定义 r = y' - f(x, y)，其中 y 是找到的解，则求解器试图在每个网格间隔上实现 norm(r / (1 + abs(f)) < tol 的标准（使用数值积分公式估计的均方根）。默认为 1e-3。

max_nodesint，可选

允许的最大网格节点数。如果超过，则算法终止。默认为 1000。

verbose{0, 1, 2}，可选

算法详细程度的级别：

• 0（默认值）：静默工作。
• 1：显示终止报告。
• 2：迭代过程中显示进展。

bc_tolfloat，可选

边界条件残差的所需绝对容差：bc 值应满足 abs(bc) < bc_tol 每个分量。默认为 tol。允许最多 10 次迭代以达到此容差。

返回：

具有以下字段定义的 Bunch 对象：

solPPoly

找到关于 y 的解为 scipy.interpolate.PPoly 实例，一个 C1 连续的三次样条。

pndarray 或 None，形状 (k,)

找到的参数。如果问题中不存在参数，则为 None。

xndarray，形状为 (m,)

最终网格的节点。

yndarray，形状为 (n, m)

在网格节点处的解值。

ypndarray，形状为 (n, m)

解在网格节点处的导数。

rms_residualsndarray，形状为 (m - 1,)

相对于每个网格间隔的相对残差的 RMS 值（请参阅 tol 参数的描述）。

niterint

完成迭代的次数。

statusint

算法终止的原因：

• 0: 算法收敛到期望的精度。
• 1: 超过了最大网格节点数。
• 2: 在解决匹配系统时遇到奇异雅可比矩阵。

messagestring

终止原因的口头描述。

successbool

如果算法收敛到期望的精度（status=0）则返回真。

注释

此函数实现了一个具有残差控制的 4 阶匹配算法，类似于 [1]。一个匹配系统通过一个具有仿射不变判据函数的阻尼牛顿法解决，如 [3] 所述。

注意，在 [1] 中，积分残差的定义没有通过区间长度进行归一化。因此，它们的定义与此处使用的定义相差一个 h**0.5 的乘数（其中 h 是区间长度）。

从版本 0.18.0 开始新增。

参考文献

[1] (1,2)

J. Kierzenka, L. F. Shampine, “基于残差控制和 Maltab PSE 的 BVP 求解器”，ACM 数学软件交易，第 27 卷，第 3 期，2001 年，299-316 页。

[2]

L.F. Shampine, P. H. Muir 和 H. Xu，“一个用户友好的 Fortran BVP 求解器”。

[3]

U. Ascher, R. Mattheij 和 R. Russell，“常微分方程边值问题的数值解法”。

[4]

共轭复数-黎曼方程 在维基百科上。

示例

在第一个例子中，我们解决 Bratu 的问题：

y'' + k * exp(y) = 0
y(0) = y(1) = 0 

对于 k = 1。

我们将方程重写为一个一阶系统，并实现其右手边的评估：

y1' = y2
y2' = -exp(y1) 

>>> import numpy as np
>>> def fun(x, y):
...     return np.vstack((y[1], -np.exp(y[0]))) 

实现边界条件残差的评估：

>>> def bc(ya, yb):
...     return np.array([ya[0], yb[0]]) 

定义具有 5 个节点的初始网格：

>>> x = np.linspace(0, 1, 5) 

这个问题已知有两个解。为了获得这两个解，我们对 y 使用两个不同的初始猜测，分别用下标 a 和 b 表示。

>>> y_a = np.zeros((2, x.size))
>>> y_b = np.zeros((2, x.size))
>>> y_b[0] = 3 

现在我们准备运行求解器。

>>> from scipy.integrate import solve_bvp
>>> res_a = solve_bvp(fun, bc, x, y_a)
>>> res_b = solve_bvp(fun, bc, x, y_b) 

让我们绘制这两个找到的解。我们利用解的样条形式来产生平滑的图形。

>>> x_plot = np.linspace(0, 1, 100)
>>> y_plot_a = res_a.sol(x_plot)[0]
>>> y_plot_b = res_b.sol(x_plot)[0]
>>> import matplotlib.pyplot as plt
>>> plt.plot(x_plot, y_plot_a, label='y_a')
>>> plt.plot(x_plot, y_plot_b, label='y_b')
>>> plt.legend()
>>> plt.xlabel("x")
>>> plt.ylabel("y")
>>> plt.show() 

我们看到两个解的形状相似，但在尺度上有显著差异。

在第二个例子中，我们解决一个简单的 Sturm-Liouville 问题：

y'' + k**2 * y = 0
y(0) = y(1) = 0 

已知对于 k = pi * n（其中 n 是整数），存在非平凡解 y = A * sin(k * x)。为了建立归一化常数 A = 1，我们添加一个边界条件：

y'(0) = k 

再次，我们将方程重写为一个一阶系统，并实现其右手边的评估：

y1' = y2
y2' = -k**2 * y1 

>>> def fun(x, y, p):
...     k = p[0]
...     return np.vstack((y[1], -k**2 * y[0])) 

注意，在 [1] 中，参数 p 被作为一个向量传递（在我们的情况下只有一个元素）。

实现边界条件：

>>> def bc(ya, yb, p):
...     k = p[0]
...     return np.array([ya[0], yb[0], ya[1] - k]) 

设置初始网格和 y 的猜测。我们旨在找到 k = 2 * pi 的解，为此我们设置 y 的值以近似 sin(2 * pi * x)：

>>> x = np.linspace(0, 1, 5)
>>> y = np.zeros((2, x.size))
>>> y[0, 1] = 1
>>> y[0, 3] = -1 

使用 6 作为 k 的初始猜测来运行求解器。

>>> sol = solve_bvp(fun, bc, x, y, p=[6]) 

我们看到找到的 k 大致正确：

>>> sol.p[0]
6.28329460046 

最后，绘制解以查看预期的正弦波形：

>>> x_plot = np.linspace(0, 1, 100)
>>> y_plot = sol.sol(x_plot)[0]
>>> plt.plot(x_plot, y_plot)
>>> plt.xlabel("x")
>>> plt.ylabel("y")
>>> plt.show() 

插值（scipy.interpolate）

原文链接：docs.scipy.org/doc/scipy-1.12.0/reference/interpolate.html

用于插值的子包对象。

如下所列，这个子包含有样条函数和类、1-D 和多维（单变量和多变量）插值类、Lagrange 和 Taylor 多项式插值器，以及 FITPACK 和 DFITPACK 函数的封装。

单变量插值

interp1d(x, y[, kind, axis, copy, …])	对一个 1-D 函数进行插值
BarycentricInterpolator(xi[, yi, axis, wi, …])	一组点的插值多项式
KroghInterpolator(xi, yi[, axis])	一组点的插值多项式
barycentric_interpolate(xi, yi, x[, axis, der])	多项式插值的便捷函数
krogh_interpolate(xi, yi, x[, der, axis])	多项式插值的便捷函数
pchip_interpolate(xi, yi, x[, der, axis])	pchip 插值的便捷函数
CubicHermiteSpline(x, y, dydx[, axis, …])	分段立方插值器，匹配值和一阶导数
PchipInterpolator(x, y[, axis, extrapolate])	PCHIP 1-D 单调立方插值器
Akima1DInterpolator(x, y[, axis])	Akima 插值器
CubicSpline(x, y[, axis, bc_type, extrapolate])	立方样条数据插值器
PPoly(c, x[, extrapolate, axis])	由系数和断点表示的分段多项式
BPoly(c, x[, extrapolate, axis])	根据系数和分断点的分段多项式。

多变量插值

非结构化数据：

griddata(points, values, xi[, method, …])	插值非结构化 D-D 数据。
LinearNDInterpolator(points, values[, …])	N > 1 维的分段线性插值器。
NearestNDInterpolator(x, y[, rescale, …])	最近邻插值器 NearestNDInterpolator(x, y)。
CloughTocher2DInterpolator(points, values[, …])	CloughTocher2DInterpolator(points, values, tol=1e-6)。
RBFInterpolator(y, d[, neighbors, …])	N 维径向基函数（RBF）插值。
Rbf(*args, **kwargs)	从 N-D 离散数据到 M-D 域的径向基函数插值的类。
interp2d(x, y, z[, kind, copy, …])

自版本 1.10.0 起不推荐使用。

|

对于网格数据：

interpn(points, values, xi[, method, …])	在正则或直角网格上进行多维插值。
RegularGridInterpolator(points, values[, …])	在任意维度的正则或直角网格上的插值器。
RectBivariateSpline(x, y, z[, bbox, kx, ky, s])	矩形网格上的双变量样条插值。

See also

scipy.ndimage.map_coordinates

张量积多项式：

NdPPoly(c, x[, extrapolate])	分段张量积多项式
NdBSpline(t, c, k, *[, extrapolate])	张量积样条对象。

1-D Splines

BSpline(t, c, k[, extrapolate, axis])	B 样条基础上的单变量样条。
make_interp_spline(x, y[, k, t, bc_type, …])	计算（系数的）插值 B 样条。
make_lsq_spline(x, y, t[, k, w, axis, …])	计算（系数的）基于 LSQ（最小二乘）拟合的 B 样条。
make_smoothing_spline(x, y[, w, lam])	使用 lam 控制曲线平滑度和数据接近度之间的权衡，计算（系数的）平滑立方样条函数。

Functional interface to FITPACK routines:

splrep(x, y[, w, xb, xe, k, task, s, t, …])	查找 1-D 曲线的 B 样条表示。
splprep(x[, w, u, ub, ue, k, task, s, t, …])	查找 N-D 曲线的 B 样条表示。
splev(x, tck[, der, ext])	计算 B 样条或其导数的值。
splint(a, b, tck[, full_output])	计算 B 样条在给定两点之间的定积分。
sproot(tck[, mest])	查找立方 B 样条的根。
spalde(x, tck)	计算 B 样条的所有导数。
splder(tck[, n])	计算给定样条导数的样条表示。
splantider(tck[, n])	计算给定样条的反导数（积分）样条。
insert(x, tck[, m, per])	在 B 样条中插入结点。

Object-oriented FITPACK interface:

UnivariateSpline(x, y[, w, bbox, k, s, ext, …])	给定数据点的一维平滑样条拟合。
InterpolatedUnivariateSpline(x, y[, w, …])	给定数据点的一维插值样条。
LSQUnivariateSpline(x, y, t[, w, bbox, k, …])	具有显式内结的一维样条。

2-D Splines

网格数据：

RectBivariateSpline(x, y, z[, bbox, kx, ky, s])	矩形网格上的二元样条逼近。
RectSphereBivariateSpline(u, v, r[, s, …])	球面上矩形网格的二元样条逼近。

无结构数据：

BivariateSpline()	二元样条基类。
SmoothBivariateSpline(x, y, z[, w, bbox, …])	平滑二元样条逼近。
SmoothSphereBivariateSpline(theta, phi, r[, …])	球坐标中平滑二元样条逼近。
LSQBivariateSpline(x, y, z, tx, ty[, w, …])	加权最小二乘二元样条逼近。
LSQSphereBivariateSpline(theta, phi, r, tt, tp)	球坐标中加权最小二乘二元样条逼近。

FITPACK 函数的底层接口：

bisplrep(x, y, z[, w, xb, xe, yb, ye, kx, …])	表面的二元 B 样条表示。
bisplev(x, y, tck[, dx, dy])	计算二元 B 样条及其导数。

附加工具

lagrange(x, w)	返回拉格朗日插值多项式。
approximate_taylor_polynomial(f, x, degree, …)	通过多项式拟合估算 f 在 x 处的 Taylor 多项式。
pade(an, m[, n])	返回多项式的 Pade 近似作为两个多项式的比值。

另请参阅

scipy.ndimage.map_coordinates, scipy.ndimage.spline_filter, scipy.signal.resample, scipy.signal.bspline, scipy.signal.gauss_spline, scipy.signal.qspline1d, scipy.signal.cspline1d, scipy.signal.qspline1d_eval, scipy.signal.cspline1d_eval, scipy.signal.qspline2d, scipy.signal.cspline2d.

pchip 是 PchipInterpolator 的别名，用于向后兼容性（新代码中不应使用）。

scipy.interpolate.interp1d

原文：docs.scipy.org/doc/scipy-1.12.0/reference/generated/scipy.interpolate.interp1d.html#scipy.interpolate.interp1d

class scipy.interpolate.interp1d(x, y, kind='linear', axis=-1, copy=True, bounds_error=None, fill_value=nan, assume_sorted=False)

插值 1-D 函数。

遗留版本

此类被视为遗留版本，将不再接收更新。这也可能意味着在未来的 SciPy 版本中将其移除。有关预期替代 interp1d 的指南，请参阅 1-D interpolation。

x 和 y 是用来近似某个函数 f 的值的数组：y = f(x)。此类返回一个函数，其调用方法使用插值来找到新点的值。

参数：

x(npoints, ) 类似数组

一个包含实数值的 1-D 数组。

y(…, npoints, …) 类似数组

一个包含实数值的 N-D 数组。沿插值轴的 y 的长度必须等于 x 的长度。使用 axis 参数选择正确的轴。与其他插值器不同，默认插值轴是 y 的最后一个轴。

kindstr 或 int，可选

指定插值类型的字符串或指定要使用的样条插值器的顺序的整数。字符串必须是以下之一：‘linear’、‘nearest’、‘nearest-up’、‘zero’、‘slinear’、‘quadratic’、‘cubic’、‘previous’ 或 ‘next’。‘zero’、‘slinear’、‘quadratic’ 和 ‘cubic’ 分别指零阶、一阶、二阶或三阶样条插值；‘previous’ 和 ‘next’ 分别返回点的前一个或后一个值；‘nearest-up’ 和 ‘nearest’ 在插值半整数（例如 0.5、1.5）时有所不同，‘nearest-up’ 向上取整，而 ‘nearest’ 向下取整。默认为 ‘linear’。

axisint，可选

y 数组中对应于 x 坐标值的轴。与其他插值器不同，默认为 axis=-1。

copybool，可选

如果为 True，则该类会对 x 和 y 进行内部复制。如果为 False，则使用对 x 和 y 的引用。默认情况下为复制。

bounds_errorbool，可选

如果为 True，在尝试在 x 范围之外的值（需要外推）进行插值时会引发 ValueError。如果为 False，则将超出范围的值分配给 fill_value。默认情况下会引发错误，除非 fill_value="extrapolate"。

fill_value类似数组或（类似数组，类似数组）或“extrapolate”，可选

• 如果是 ndarray（或 float），则此值将用于填充数据范围外请求的点。如果未提供，则默认为 NaN。类似数组必须正确传播到非插值轴的维度。

• 如果是两个元素的元组，则第一个元素用作x_new < x[0]的填充值，第二个元素用于x_new > x[-1]。任何不是两个元素的元组（例如列表或 ndarray，无论形状如何）都被视为单个类似数组的参数，用于below, above = fill_value, fill_value。使用两个元素的元组或 ndarray 需要bounds_error=False。

  版本 0.17.0 中的新功能。

• 如果为“extrapolate”，则数据范围外的点将被外推。

  版本 0.17.0 中的新功能。

assume_sortedbool，可选

如果为 False，则x的值可以按任意顺序排列，并且它们首先被排序。如果为 True，则x必须是单调递增值的数组。

另请参阅

splrep，splev

基于 FITPACK 的样条插值/平滑。

UnivariateSpline

FITPACK 例程的面向对象封装。

interp2d

二维插值

注意事项

在输入值中存在 NaN 时调用interp1d将导致未定义的行为。

输入值x和y必须可转换为float值，如int或float。

如果x中的值不唯一，则结果行为未定义，并且取决于kind的选择，即更改kind会改变重复项的行为。

示例

>>> import numpy as np
>>> import matplotlib.pyplot as plt
>>> from scipy import interpolate
>>> x = np.arange(0, 10)
>>> y = np.exp(-x/3.0)
>>> f = interpolate.interp1d(x, y) 

>>> xnew = np.arange(0, 9, 0.1)
>>> ynew = f(xnew)   # use interpolation function returned by `interp1d`
>>> plt.plot(x, y, 'o', xnew, ynew, '-')
>>> plt.show() 

属性：

fill_value

填充值。

方法

__call__(x)	评估插值

scipy.interpolate.BarycentricInterpolator

原文：docs.scipy.org/doc/scipy-1.12.0/reference/generated/scipy.interpolate.BarycentricInterpolator.html#scipy.interpolate.BarycentricInterpolator

class scipy.interpolate.BarycentricInterpolator(xi, yi=None, axis=0, *, wi=None, random_state=None)

一组点的插值多项式。

构造通过给定点集的多项式。允许评估多项式及其所有导数，有效更改要插值的 y 值，并通过添加更多的 x 和 y 值进行更新。

出于数值稳定性的原因，此函数不计算多项式的系数。

在函数评估之前需要提供 yi 的值，但前处理不依赖于它们，因此可以快速更新。

参数：

xiarray_like，形状为 (npoints, )

一维数组，多项式应通过的点的 x 坐标

yiarray_like，形状为 (…, npoints, …)，可选

y 坐标的 N-D 数组，多项式应通过这些点。如果为 None，则稍后通过 set_y 方法提供 y 值。沿插值轴的 yi 的长度必须等于 xi 的长度。使用 axis 参数选择正确的轴。

axisint，可选

yi 数组中对应于 x 坐标值的轴。默认为 axis=0。

wiarray_like，可选

所选插值点 xi 的重心权重。如果缺少或为 None，则从 xi 计算权重（默认）。这允许在使用相同节点 xi 计算多个插值时重复使用权重 wi，而无需重新计算。

random_state{None, int, numpy.random.Generator，numpy.random.RandomState}，可选

如果 seed 为 None（或 np.random），则使用 numpy.random.RandomState 单例。如果 seed 是整数，则使用新的 RandomState 实例，并使用 seed 进行种子化。如果 seed 已经是 Generator 或 RandomState 实例，则使用该实例。

注意事项

此类使用“重心插值”方法，将问题视为有理函数插值的特例。这种算法在数值上非常稳定，但即使在精确计算的世界中，除非选择 x 坐标非常仔细 - Chebyshev 零点（例如，cos(i*pi/n)）是一个很好的选择 - 多项式插值本身也是一个非常病态的过程，这是由于 Runge 现象。

基于 Berrut 和 Trefethen 2004 年的“Barycentric Lagrange Interpolation”。

示例

要在区间 ((0, \frac{\pi}{2})) 内使用六个随机分布的节点生成一个近似于函数 (\sin x) 及其前四阶导数的五次重心插值多项式：

>>> import numpy as np
>>> import matplotlib.pyplot as plt
>>> from scipy.interpolate import BarycentricInterpolator
>>> rng = np.random.default_rng()
>>> xi = rng.random(6) * np.pi/2
>>> f, f_d1, f_d2, f_d3, f_d4 = np.sin, np.cos, lambda x: -np.sin(x), lambda x: -np.cos(x), np.sin
>>> P = BarycentricInterpolator(xi, f(xi), random_state=rng)
>>> fig, axs = plt.subplots(5, 1, sharex=True, layout='constrained', figsize=(7,10))
>>> x = np.linspace(0, np.pi, 100)
>>> axs[0].plot(x, P(x), 'r:', x, f(x), 'k--', xi, f(xi), 'xk')
>>> axs[1].plot(x, P.derivative(x), 'r:', x, f_d1(x), 'k--', xi, f_d1(xi), 'xk')
>>> axs[2].plot(x, P.derivative(x, 2), 'r:', x, f_d2(x), 'k--', xi, f_d2(xi), 'xk')
>>> axs[3].plot(x, P.derivative(x, 3), 'r:', x, f_d3(x), 'k--', xi, f_d3(xi), 'xk')
>>> axs[4].plot(x, P.derivative(x, 4), 'r:', x, f_d4(x), 'k--', xi, f_d4(xi), 'xk')
>>> axs[0].set_xlim(0, np.pi)
>>> axs[4].set_xlabel(r"$x$")
>>> axs[4].set_xticks([i * np.pi / 4 for i in range(5)],
...                   ["0", r"$\frac{\pi}{4}$", r"$\frac{\pi}{2}$", r"$\frac{3\pi}{4}$", r"$\pi$"])
>>> axs[0].set_ylabel("$f(x)$")
>>> axs[1].set_ylabel("$f'(x)$")
>>> axs[2].set_ylabel("$f''(x)$")
>>> axs[3].set_ylabel("$f^{(3)}(x)$")
>>> axs[4].set_ylabel("$f^{(4)}(x)$")
>>> labels = ['Interpolation nodes', 'True function $f$', 'Barycentric interpolation']
>>> axs[0].legend(axs[0].get_lines()[::-1], labels, bbox_to_anchor=(0., 1.02, 1., .102),
...               loc='lower left', ncols=3, mode="expand", borderaxespad=0., frameon=False)
>>> plt.show() 

属性:

dtype

方法

__call__(x)	在点 x 处评估插值多项式
add_xi(xi[, yi])	将更多 x 值添加到待插值的集合中
derivative(x[, der])	在点 x 处评估多项式的单个导数
derivatives(x[, der])	在点 x 处评估多项式的多个导数
set_yi(yi[, axis])	更新待插值的 y 值

scipy.interpolate.KroghInterpolator

原文链接：docs.scipy.org/doc/scipy-1.12.0/reference/generated/scipy.interpolate.KroghInterpolator.html#scipy.interpolate.KroghInterpolator

class scipy.interpolate.KroghInterpolator(xi, yi, axis=0)

一组点的插值多项式。

该多项式通过所有配对的(xi, yi)。还可以指定每个点xi处的多个导数；通过重复值xi并按顺序指定导数值yi来完成。

允许评估多项式及其所有导数。出于数值稳定性的原因，此函数不计算多项式的系数，但可以通过评估所有导数来获得它们。

参数：

xi类数组，形状（npoints，）

已知的 x 坐标。必须按升序排列。

yi类数组，形状（…，npoints，…）

已知的 y 坐标。当 xi 连续出现两次或更多时，对应的 yi 表示导数值。沿插值轴的yi的长度必须等于xi的长度。使用axis参数选择正确的轴。

axis整型，可选

在yi数组中对应于 x 坐标值的轴。默认为axis=0。

注意事项

请注意，这里实现的算法不一定是已知的最稳定的。此外，即使在精确计算的世界中，除非选择的 x 坐标非常谨慎 - Chebyshev 零点（例如，cos(i*pi/n)）是一个很好的选择 - 多项式插值本身也是一个非常病态的过程，因为 Runge 现象。一般来说，即使选择了良好的 x 值，在本代码中，大于约 30 的度数会导致数值不稳定性问题。

基于[1]。

参考资料

[1]

Krogh，《多项式插值和数值微分的高效算法》，1970 年。

示例

要生成一个在 0 和 1 处为零且在 0 处导数为 2 的多项式，请调用

>>> from scipy.interpolate import KroghInterpolator
>>> KroghInterpolator([0,0,1],[0,2,0]) 

这构造了二次多项式(2x²-2x)。在xi数组中通过重复的零指示导数条件；对应的 yi 值为 0，函数值为 2，导数值为 2。

举个例子，对于给定的xi、yi和每个点的导数ypi，可以构建适当的数组如下：

>>> import numpy as np
>>> rng = np.random.default_rng()
>>> xi = np.linspace(0, 1, 5)
>>> yi, ypi = rng.random((2, 5))
>>> xi_k, yi_k = np.repeat(xi, 2), np.ravel(np.dstack((yi,ypi)))
>>> KroghInterpolator(xi_k, yi_k) 

要生成一个向量值多项式，请为yi提供一个更高维度的数组：

>>> KroghInterpolator([0,1],[[2,3],[4,5]]) 

这构造了一个线性多项式，在 0 处给出(2,3)，在 1 处给出(4,5)。

属性：

dtype

方法

__call__(x)	评估插值
derivative(x[, der])	在点 x 处评估单个多项式导数。
derivatives(x[, der])	在点 x 处评估多个多项式导数。

scipy.interpolate.barycentric_interpolate

原文：docs.scipy.org/doc/scipy-1.12.0/reference/generated/scipy.interpolate.barycentric_interpolate.html#scipy.interpolate.barycentric_interpolate

scipy.interpolate.barycentric_interpolate(xi, yi, x, axis=0, *, der=0)

多项式插值的便捷函数。

构造一个通过给定点集的多项式，然后评估多项式。由于数值稳定性的原因，此函数不计算多项式的系数。

此函数使用“重心插值”方法，将问题视为有理函数插值的特殊情况。这种算法在数值上非常稳定，但即使在精确计算的世界中，除非非常谨慎地选择x坐标（例如，切比雪夫零点（例如，cos(i*pi/n)）是一个很好的选择），多项式插值本身也是一个由于 Runge 现象而非常病态的过程。

参数：

xiarray_like

多项式应通过的点的 x 坐标的 1-D 数组

yiarray_like

多项式应通过的点的 y 坐标。

x标量 或 array_like

评估插值函数的点或点。

derint 或 列表 或 None, 可选

要评估的导数数量，或者对所有可能非零导数评估（即与点数相等的数量），或者要评估的导数列表。这个数字包括函数值作为“0th”导数。

axisint, 可选

与yi数组中的轴相对应的 x 坐标值。

返回：

y标量 或 array_like

插值值。形状由用x替换原始数组中的插值轴决定。

另请参阅

BarycentricInterpolator

重心插值器

注解

插值权重的构造是一个相对较慢的过程。如果您希望多次使用相同的 xi 调用此函数（但可能会变化 yi 或 x），您应该使用类BarycentricInterpolator。这是此函数在内部使用的内容。

示例

我们可以使用重心插值法对 2D 观察数据进行插值：

>>> import numpy as np
>>> import matplotlib.pyplot as plt
>>> from scipy.interpolate import barycentric_interpolate
>>> x_observed = np.linspace(0.0, 10.0, 11)
>>> y_observed = np.sin(x_observed)
>>> x = np.linspace(min(x_observed), max(x_observed), num=100)
>>> y = barycentric_interpolate(x_observed, y_observed, x)
>>> plt.plot(x_observed, y_observed, "o", label="observation")
>>> plt.plot(x, y, label="barycentric interpolation")
>>> plt.legend()
>>> plt.show() 

scipy.interpolate.krogh_interpolate

原文：docs.scipy.org/doc/scipy-1.12.0/reference/generated/scipy.interpolate.krogh_interpolate.html#scipy.interpolate.krogh_interpolate

scipy.interpolate.krogh_interpolate(xi, yi, x, der=0, axis=0)

用于多项式插值的便捷函数。

参见 KroghInterpolator 了解更多细节。

参数：

xiarray_like

插值点（已知 x 坐标）。

yiarray_like

已知的 y 坐标，形状为 (xi.size, R)。如果 R=1，则解释为长度为 R 的向量或标量。

xarray_like

要评估导数的点或点。

derint 或 列表 或 None，可选

要评估的导数数量，或者对所有可能非零导数（即与点数相等的数字）进行评估，或者要评估的导数列表。该数字包括函数值作为第 ‘0’ 导数。

axisint，可选

yi 数组中对应于 x 坐标值的轴。

返回：

dndarray

如果插值器的值为 R-D，则返回的数组将为 N by R 的导数数量。如果 x 是标量，则将去掉中间维度；如果 yi 是标量，则将去掉最后维度。

参见

KroghInterpolator

Krogh 插值器

注意

插值多项式的构造是一个相对昂贵的过程。如果需要重复评估它，请考虑使用类 KroghInterpolator（这正是该函数使用的内容）。

示例

我们可以使用 Krogh 插值来插值 2D 观测数据：

>>> import numpy as np
>>> import matplotlib.pyplot as plt
>>> from scipy.interpolate import krogh_interpolate
>>> x_observed = np.linspace(0.0, 10.0, 11)
>>> y_observed = np.sin(x_observed)
>>> x = np.linspace(min(x_observed), max(x_observed), num=100)
>>> y = krogh_interpolate(x_observed, y_observed, x)
>>> plt.plot(x_observed, y_observed, "o", label="observation")
>>> plt.plot(x, y, label="krogh interpolation")
>>> plt.legend()
>>> plt.show() 

import numpy as np

import matplotlib.pyplot as plt
from scipy.interpolate import barycentric_interpolate
x_observed = np.linspace(0.0, 10.0, 11)
y_observed = np.sin(x_observed)
x = np.linspace(min(x_observed), max(x_observed), num=100)
y = barycentric_interpolate(x_observed, y_observed, x)
plt.plot(x_observed, y_observed, “o”, label=“observation”)
plt.plot(x, y, label=“barycentric interpolation”)
plt.legend()
plt.show()

[外链图片转存中...(img-LmlfmCXa-1719633858086)]


# `scipy.interpolate.krogh_interpolate`

> 原文：[`docs.scipy.org/doc/scipy-1.12.0/reference/generated/scipy.interpolate.krogh_interpolate.html#scipy.interpolate.krogh_interpolate`](https://docs.scipy.org/doc/scipy-1.12.0/reference/generated/scipy.interpolate.krogh_interpolate.html#scipy.interpolate.krogh_interpolate)

```py
scipy.interpolate.krogh_interpolate(xi, yi, x, der=0, axis=0)

用于多项式插值的便捷函数。

参见 KroghInterpolator 了解更多细节。

参数：

xiarray_like

插值点（已知 x 坐标）。

yiarray_like

已知的 y 坐标，形状为 (xi.size, R)。如果 R=1，则解释为长度为 R 的向量或标量。

xarray_like

要评估导数的点或点。

derint 或 列表 或 None，可选

要评估的导数数量，或者对所有可能非零导数（即与点数相等的数字）进行评估，或者要评估的导数列表。该数字包括函数值作为第 ‘0’ 导数。

axisint，可选

yi 数组中对应于 x 坐标值的轴。

返回：

dndarray

如果插值器的值为 R-D，则返回的数组将为 N by R 的导数数量。如果 x 是标量，则将去掉中间维度；如果 yi 是标量，则将去掉最后维度。

参见

KroghInterpolator

Krogh 插值器

注意

插值多项式的构造是一个相对昂贵的过程。如果需要重复评估它，请考虑使用类 KroghInterpolator（这正是该函数使用的内容）。

示例

我们可以使用 Krogh 插值来插值 2D 观测数据：

>>> import numpy as np
>>> import matplotlib.pyplot as plt
>>> from scipy.interpolate import krogh_interpolate
>>> x_observed = np.linspace(0.0, 10.0, 11)
>>> y_observed = np.sin(x_observed)
>>> x = np.linspace(min(x_observed), max(x_observed), num=100)
>>> y = krogh_interpolate(x_observed, y_observed, x)
>>> plt.plot(x_observed, y_observed, "o", label="observation")
>>> plt.plot(x, y, label="krogh interpolation")
>>> plt.legend()
>>> plt.show() 

[外链图片转存中…(img-hZsYoqf3-1719633858086)]