Issue 1:
I am having an issue with what seems to be a very infrequent infinite loop in RungeKutta._step_impl in scipy/integrate/_ivp/rk.py. Running a bunch of experiments in parallel on large number of cores, I maybe see it once a week, and not in any easily reproducible fashion, as far as I can tell.

I caught it after, as far as I can tell, scipy.integrate was running for a couple of days without completing. I am guessing RungeKutta._step_impl is the one not completing, but perhaps it's higher in the call stack.

I examined the running process using gdb with python extensions. At the point where rk_step is called inside the while loop, the following are some of the relevant values (more details are in the gdb session log below):

self.n == 6
self.t == t == 0
self.y == y == numpy.ndarray([2.19183852e-33, 0, 1.40129846e-45, 0, 1.39532281e-34, 0])
self.h_abs == -NaN
h_abs == NaN
self.max_step == max_step == Inf
min_step == 4.9406564584124654e-323
self.direction == 1
order == 4
h == t_new == -NaN
self.t_bound = 0.02

I am running:

Ubuntu 16.04

python 3.5.2

numpy 1.15.4

scipy 1.3.3

More complete gdb output:

(gdb) py-bt
Traceback (most recent call first):
[...]
File "[...]/envs/model/model.py", line 215, in _tire_dynamics_rear
* (np.tan(alpha)/(1+K))**2 )
File "[...]/envs/model/model.py", line 83, in _dynamics
F_xr, F_yr = self._tire_dynamics_rear(v_x, cmd_vx, alpha_r)
File "[...]/envs/model/model.py", line 55, in <lambda>
fun=(lambda t, X: self._dynamics(X, t, U)),
File "/usr/local/lib/python3.5/dist-packages/scipy/integrate/_ivp/base.py", line 21, in fun_wrapped
return np.asarray(fun(t, y), dtype=dtype)
File "/usr/local/lib/python3.5/dist-packages/scipy/integrate/_ivp/base.py", line 139, in fun
return self.fun_single(t, y)
File "/usr/local/lib/python3.5/dist-packages/scipy/integrate/_ivp/rk.py", line 70, in rk_step
K[s + 1] = fun(t + c * h, y + dy)
File "/usr/local/lib/python3.5/dist-packages/scipy/integrate/_ivp/rk.py", line 143, in _step_impl
self.B, self.C, self.E, self.K)
File "/usr/local/lib/python3.5/dist-packages/scipy/integrate/_ivp/base.py", line 182, in step
success, message = self._step_impl()
File "/usr/local/lib/python3.5/dist-packages/scipy/integrate/_ivp/ivp.py", line 502, in solve_ivp
message = solver.step()
File "[...]/envs/model/model.py", line 56, in state_transition
t_span=t, y0=X, atol=1e-5)
[...]
(gdb) py-up
[...]
(gdb) py-up
#24 Frame 0x61ffa18, for file /usr/local/lib/python3.5/dist-packages/scipy/integrate/_ivp/rk.py, line 143, in _step_impl (self=<RK45(t=<float at remote 0x7f80dcc71510>, direction=<numpy.float64 at remote 0x7f80dccc5ee8>, K=<numpy.ndarray at remote 0x7f80dcc53d00>, vectorized=False, _fun=<function at remote 0x7f80dcc51f28>, y_old=None, status='running', f=<numpy.ndarray at remote 0x7f80dcc53a80>, y=<numpy.ndarray at remote 0x7f80dcc53350>, fun_vectorized=<function at remote 0x7f80dcc51950>, t_old=None, fun=<function at remote 0x7f80dcc51d90>, t_bound=<float at remote 0x7f80dcc71090>, njev=0, fun_single=<function at remote 0x7f80dcc51f28>, nfev=2149842846, nlu=0, atol=<numpy.ndarray at remote 0x7f80dcc537b0>, h_abs=<numpy.float64 at remote 0x7f80dccdb2d0>, rtol=<float at remote 0x7f809c0ca6a8>, n=6, max_step=<float at remote 0x7f8184eb3420>) at remote 0x7f80a0c9bb00>, t=<float at remote 0x7f80dcc71510>, y=<numpy.ndarray at remote 0x7f80dcc53350>, max_step=<float at remote 0x7f8184eb3420>, rtol=<float at remote 0x7f...(truncated)
self.B, self.C, self.E, self.K)
(gdb) py-print self
local 'self' = <RK45(t=<float at remote 0x7f80dcc71510>, direction=<numpy.float64 at remote 0x7f80dccc5ee8>, K=<numpy.ndarray at remote 0x7f80dcc53d00>, vectorized=False, _fun=<function at remote 0x7f80dcc51f28>, y_old=None, status='running', f=<numpy.ndarray at remote 0x7f80dcc53a80>, y=<numpy.ndarray at remote 0x7f80dcc53350>, fun_vectorized=<function at remote 0x7f80dcc51950>, t_old=None, fun=<function at remote 0x7f80dcc51d90>, t_bound=<float at remote 0x7f80dcc71090>, njev=0, fun_single=<function at remote 0x7f80dcc51f28>, nfev=2149842846, nlu=0, atol=<numpy.ndarray at remote 0x7f80dcc537b0>, h_abs=<numpy.float64 at remote 0x7f80dccdb2d0>, rtol=<float at remote 0x7f809c0ca6a8>, n=6, max_step=<float at remote 0x7f8184eb3420>) at remote 0x7f80a0c9bb00>
(gdb) py-print t
local 't' = <float at remote 0x7f80dcc71510>
(gdb) p * (PyFloatObject*) 0x7f80dcc71510
$36 = {ob_base = {ob_refcnt = 7, ob_type = 0xa3c2c0 <PyFloat_Type>}, ob_fval = 0}
(gdb) py-print y
local 'y' = <numpy.ndarray at remote 0x7f80dcc53350>
(gdb) p *(PyArrayObject_fields*) 0x7f80dcc53350
$37 = {ob_base = {ob_refcnt = 12, ob_type = 0x7f81838f7180},
data = 0x9361740 "4D\350\372^\a\005\300\320:\240v<M\374?\221\361\240\266#\343\347\277\250\037\335d\323\271\376?O\210g\377\221\331\307?\332\v+1\227\315\301?\320\246\235\006", nd = 1, dimensions = 0x7397860,
strides = 0x7397868, base = <numpy.ndarray at remote 0x7f80dcc53710>, descr = 0x7f81838f7ec0, flags = 1283,
weakreflist = 0x0}
(gdb) p ((PyArrayObject_fields*) 0x7f80dcc53350)->dimensions
$38 = (npy_intp *) 0x7397860
(gdb)
$39 = (npy_intp *) 0x7397860
(gdb) p *((PyArrayObject_fields*) 0x7f80dcc53350)->dimensions
$40 = 6
(gdb) p (float[6])(((PyArrayObject_fields*) 0x7f80dcc53350)->data)
$41 = {2.19183852e-33, 0, 1.40129846e-45, 0, 1.39532281e-34, 0}
(gdb) py-print h_abs
local 'h_abs' = <numpy.float64 at remote 0x7f80dccdbca8>
(gdb) p * (PyFloatObject*) 0x7f80dccdb2d0
$42 = {ob_base = {ob_refcnt = 1, ob_type = 0x7f8183909600}, ob_fval = -nan(0x8000000000000)}
(gdb) p * (PyFloatObject*) 0x7f80dccdbca8
$43 = {ob_base = {ob_refcnt = 1, ob_type = 0x7f8183909600}, ob_fval = nan(0x8000000000000)}
(gdb) py-print max_step
local 'max_step' = <float at remote 0x7f8184eb3420>
(gdb) p * (PyFloatObject*) 0x7f8184eb3420
$44 = {ob_base = {ob_refcnt = 246, ob_type = 0xa3c2c0 <PyFloat_Type>}, ob_fval = inf}
(gdb) py-print min_step
local 'min_step' = <numpy.float64 at remote 0x7f80dcd0f360>
(gdb) p * (PyFloatObject*) 0x7f80dcd0f360
$45 = {ob_base = {ob_refcnt = 1, ob_type = 0x7f8183909600}, ob_fval = 4.9406564584124654e-323}
(gdb) p * (PyFloatObject*) 0x7f80dccc5ee8
$46 = {ob_base = {ob_refcnt = 1, ob_type = 0x7f8183909600}, ob_fval = 1}
(gdb) py-print order
local 'order' = 4
(gdb) py-print step_accepted
local 'step_accepted' = False
(gdb) py-print h
local 'h' = <numpy.float64 at remote 0x7f80dccd4270>
(gdb) p * (PyFloatObject*) 0x7f80dccd4270
$47 = {ob_base = {ob_refcnt = 3, ob_type = 0x7f8183909600}, ob_fval = -nan(0x8000000000000)}
(gdb) p * (PyFloatObject*) 0x7f80dcc71090
$48 = {ob_base = {ob_refcnt = 2, ob_type = 0xa3c2c0 <PyFloat_Type>}, ob_fval = 0.02}
(gdb) py-print t_new
local 't_new' = <numpy.float64 at remote 0x7f80dcdacf18>
(gdb) p * (PyFloatObject*) 0x7f80dcdacf18
$49 = {ob_base = {ob_refcnt = 1, ob_type = 0x7f8183909600}, ob_fval = -nan(0x8000000000000)}

I am not able to share all the code, but the model.py file is essentially the same as https://github.com/r-pad/aa_simulation/blob/990f548047/envs/model/model.py

Issue 2:
Describe your issue.

scipy.integrate._ivp.common.select_initial_step() returns numpy.nan if max(abs(y0)) == numpy.inf

This causes programs to lock, when calling solve_ivp() with such an initial value.

Probably, select_initial_step() (or its callers) should raise an Error in that situation.

Reproducing Code Example

Call `solve_ivp()` with at least one component of `y0` being infinite (probably a NaN component would cause the same problem).

Error message

None, program goes into infinite loop, when calling `solve_ivp()`

SciPy/NumPy/Python version information

1.9.3 1.23.5 sys.version_info(major=3, minor=10, micro=6, releaselevel='final', serial=0)