matplotlib / scipy在计算和绘制多个数据集（ODE的相位肖像）时会产生奇怪的结果

Question

让我试着解决我遇到的问题： 给定是一个常微分方程dx / dt = f（x，t），其相位图将通过绘制初始值数组的相位轨迹来绘制，具有以下代码：

import matplotlib.pyplot as plt
import numpy as np
from scipy.integrate import odeint

def phase_trajectories(func, inits, xbound, ybound):
    """func: RHS of ODE to be integrated by odeint
    inits: list or array of initial-value-tuples
    xbound, ybound: Boundaries for plotting
    """
    t = np.linspace(0, 10, 100)  # solution is calculated for these values
    axis = plt.gca()
    def f_neg(x, t):
        return -func(x, t)
    for i in inits:
        sol_fwd = odeint(func, i, t)   # calc solution forward
        sol_bwd = odeint(f_neg, i, t)  # and backward in time
        sol = np.vstack((np.flipud(sol_bwd), sol_fwd))  # put both together
        sol_x = sol[:,0]
        sol_y = sol[:,1]
        sol_x_masked = np.ma.masked_outside(sol_x, *xbound)  # masking the data because
        sol_y_masked = np.ma.masked_outside(sol_y, *ybound)  # of blow-up of the solution
        axis.plot(sol_x_masked, sol_y_masked, c='b')
    plt.plot()

作为一个具体的例子，看看dx / dt = x（t）²，初始值为x（t0）= x0：

# inits contains tuples (t0, x0)
inits = np.array([(i, j) for i in range(-2, 3) for j in (-1, 1)])

def func(x, t):
    # receives x as array with shape (n,) from odeint
    # returns [dt/dt, dx/dt]
    return np.array([1, x[1]**2])

现在调用phase_trajectories(func, inits, (-2, 2), (-2, 2))绘制预期的轨迹，但也添加了每次调用函数时似乎不同的不需要的噪声：example with noise another example

为了得到这个的底部，我调用phase_trajectories不是为了整个inits-array，而是为每个初始值对分别使用生成器手动调用它的next方法直到耗尽：

def generator_pt():
    for i in inits:
        yield phase_trajectories(func, (i,), (-2, 2), (-2, 2))

g = generator_pt()
next(g)

然而，这有时（奇怪的不是每次都）给我预期的结果： It's supposed to look like that

我知道给定ODE的解决方案会爆炸，这意味着它只存在于特定的间隔上并且在接近其极限时会发散，这就是为什么我屏蔽了保存数据的数组进行绘图的原因。 尽管如此，对我来说仍然很神秘，为什么每次计算和绘制数据时不需要的噪声都不同。

基于此处陈述的观察结果，我认为这种不当行为的原因是要在matplotlib中寻找，或者更确切地说是在我这方面的一些可怕的误用中。我已经尝试过使用不同的matplotlib-backends，并在不同的环境中执行代码（Python 2.7和3.5），并更新了必要的包装，但到目前为止，我的结果还没有结果，我已经达到了我只是抓住吸管的地步。

也许有人可以给我一些暗示或提供一些解释这些结果的见解，或者至少帮助我理解为什么事情会以这种方式发生。

Stackoverflow，你是我唯一的希望......

Answer 1

似乎问题在于scipys odeint方法，这是一个较老的实现，以数字方式解决ODE。从scipy v1.0.0开始，可以使用solve_ivp的新实现，不会出现干扰噪声。 以供参考：

import matplotlib.pyplot as plt
import numpy as np
from scipy.integrate import solve_ivp

def phase_trajectories(func, inits, xbound, ybound):
    axis = plt.gca()
    t = (0, 10)
    def f_neg(t, x):
        return -func(t, x)
    for f in (func, f_neg):
        for i in inits:
            x, y = solve_ivp(f, t, i).y
            x_ma = np.ma.masked_outside(x, *xbound)
            y_ma = np.ma.masked_outside(y, *ybound)
            axis.plot(x_ma, y_ma, c='b')
    plt.plot()

函数的结果和初始值looks as aspected。

然而，为什么odeint显示其奇怪行为的原因仍然没有得到答复，并且可以在底层的FORTRAN实现中找到。