встречается недопустимое значение в double_scalars def f (t, x): вернуть np.power (x, 2) - x

Question

В последнее время я много занимался построением графиков полей уклонов и решений ODE, и я решил попробовать себя в создании небольшой функции, которая автоматически отображает графики с наложением векторных полей.

Эта функция принимает набор начальных условий и строит множество решений. Это работает довольно хорошо, но для некоторых начальных значений я получаю ошибку в заголовке:

invalid value encountered in double_scalars
  def f(t, x): return np.power(x, 2) - x

Вот код функции:

def grapher(fn, t_0, t_n, dt, y_0):
"""
Takes a first order ODE and solves it for initial conditions
provided by y_0

:param fn: y' = f(t,y) 
:param t_0: start time
:param t_n: end time
:param dt:  step size
:param y_0: iterable containing initial conditions
:return: 
"""
    t = np.arange(t_0, t_n, dt)
    y_min = .0
    y_max = .0

    for iv in np.asarray(y_0):
        soln = rk4(dt, t, fn, iv)
        plt.plot(t, soln, '-r')
        if y_min > np.min(soln):
            y_min = np.min(soln)
        if y_max < np.max(soln):
            y_max = np.max(soln)

    x = np.linspace(t_0, t_n + dt, 11)
    y = np.linspace(y_min, y_max, 11)

    X, Y = np.meshgrid(x, y)

    theta = np.arctan(f(X, Y))

    U = np.cos(theta)
    V = np.sin(theta)

    plt.quiver(X, Y, U, V, angles='xy')
    plt.xlim((t_0, t_n - dt))
    plt.ylim((y_min - .1*y_min, y_max + .1*y_max))
    plt.show()

А вот приложение, котороетерпит неудачу:

def f(t, x): return x**2 - x

grapher(f,0,4,0.1, (-0.9, 0.9, 1.1))

Создает этот график, в котором отсутствует решение, связанное с начальным условием 1.1:

Однако, если я выберу значение, меньшее или равное 1, я получу правильный график:

Я не вижу возможности делить на ноль, поэтому я 'Я немного смущен. Кроме того, качественные характеристики ODE не отображаются полностью, если я не могу выбрать начальное условие выше 1.

Я хотел бы также отметить, что, когда у меня не было функции для автоматизации этого процессаопределенная функция f (x) = x ^ 2 - x не доставила мне никаких проблем. Любая подсказка, почему это может быть?

Если это поможет, вот алгоритм rk4, который я написал в другом модуле:

def rk4(dt, t, field, y_0):
    """
    :param dt: float - the timestep
    :param t: array - the time mesh
    :param field: method - the vector field y' = f(t, y)
    :param y_0: array - contains initial conditions
    :return: ndarray - solution
    """

    # Initialize solution matrix. Each row is the solution to the system
    # for a given time step. Each column is the full solution for a single
    # equation.
    y = np.asarray(len(t) * [y_0])

    for i in np.arange(len(t) - 1):
        k1 = dt * field(t[i], y[i])
        k2 = dt * field(t[i] + 0.5 * dt, y[i] + 0.5 * k1)
        k3 = dt * field(t[i] + 0.5 * dt, y[i] + 0.5 * k2)
        k4 = dt * field(t[i] + dt, y[i] + k3)
        y[i + 1] = y[i] + (k1 + 2 * k2 + 2 * k3 + k4) / 6

    return y

Answer 1

Я думаю, что в коде нет ошибки, просто решение становится слишком большим.

Если бы вы позвонили grapher с

grapher(f, 0, 4, 0.1, (-0.9, 0.9, 1.01))

, вы бы получили:

---

С:

grapher(f, 0, 4, 0.1, (-0.9, 0.9, 1.02))

и когда y_0 становится равным 1.1 значение для soln не сообщается, потому что np.pow() при обнаружении переполнения просто возвращает nan, который тогда matplotlib не знает, как построить.

---

Если вы изменили

def f(t, x):
    return x**2 - x

на:

def f(t, x):
    return x * (x - 1)

, вы бы получили (некрасивый, но "правильный") график решения для y_0 == 1.1, потому что вместо этогоиз-за переполнения по умолчанию nan теперь вы получаете inf с в качестве максимальных значений, которые, конечно, matplotlib не знает, как обрабатывать в процессе генерации осей: