Странное поведение в scipy.solve_ivp при использовании неявного метода

Question

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

from scipy.integrate import solve_ivp
import numpy as np

def f(y, t):
    return y

y0 = [1,1,1,1]
method = 'RK23'
s = solve_ivp(f, (0,1), y0, method=method, t_eval=np.linspace(0,1))

И все работает нормально.Когда я изменяю на method='BDF' или method='Radau', я получаю сообщение об ошибке:

Traceback (most recent call last):

  File "<ipython-input-222-f11c4406e92c>", line 10, in <module>
    s = solve_ivp(f, (0,1), y0, method=method, t_eval=np.linspace(0,1))

  File "C:\ProgramData\Anaconda3\lib\site-packages\scipy\integrate\_ivp\ivp.py", line 455, in solve_ivp
    solver = method(fun, t0, y0, tf, vectorized=vectorized, **options)

  File "C:\ProgramData\Anaconda3\lib\site-packages\scipy\integrate\_ivp\radau.py", line 299, in __init__
    self.jac, self.J = self._validate_jac(jac, jac_sparsity)

  File "C:\ProgramData\Anaconda3\lib\site-packages\scipy\integrate\_ivp\radau.py", line 345, in _validate_jac
    J = jac_wrapped(t0, y0, self.f)

  File "C:\ProgramData\Anaconda3\lib\site-packages\scipy\integrate\_ivp\radau.py", line 343, in jac_wrapped
    sparsity)

  File "C:\ProgramData\Anaconda3\lib\site-packages\scipy\integrate\_ivp\common.py", line 307, in num_jac
    return _dense_num_jac(fun, t, y, f, h, factor, y_scale)

  File "C:\ProgramData\Anaconda3\lib\site-packages\scipy\integrate\_ivp\common.py", line 318, in _dense_num_jac
    diff = f_new - f[:, None]

IndexError: too many indices for array

Я также получаю сообщение об ошибке с method = 'LSODA', хотя оно и отличается (т. Е. Все неявные интеграторы).Я не получаю сообщение об ошибке ни с одним из явных интеграторов.

Я пробовал это в spyder с scipy версии 1.0.0 и в google colab (scipy версии 1.1.0), с теми же результатами.

Это ошибка или мне не хватает аргумента, который мне нужен для неявных интеграторов ??