What does the letter k mean in the documentation of solve_ivp function of Scipy?

Question

Solve_ivp is an initial value problem solver function from Scipy. In a few words

scipy.integrate.solve_ivp(fun, t_span, y0, method=’RK45’, t_eval=None, dense_output=False, events=None, vectorized=False, args=None, **options)

Solve an initial value problem for a system of ODEs. This function numerically integrates a system of ordinary differential equations given an initial value.

In the solve_ivp function documentation (Scipy reference guide 1.4.1 page 695) we have the following

Parameters fun [callable] Right-hand side of the system. The calling signature is fun(t, y). Here t is a scalar, and there are two options for the ndarray y: It can either have the shape (n,); then fun must return array_like with shape (n,). Alternatively, it can have the shape (n, k); then fun must return an array_like with shape (n, k), i.e. each column corresponds to a single column in y. The choice between the two options is determined by the vectorized argument (see below). The vectorized implementation allows a faster approximation of the Jacobian by finite differences (required for stiff solvers).

Here n stands for the no of dimensions in y. What does k stand for? This might seem a very naive question for those who know the answer. But please, believe me, I really could not find it (at least not in the documentation). The great hpaulj's answer to this question seems to shed some light. But well, IMHO it still is too dark to move around.

Frazier

This is still a good question, if you've found an answer it's perfectly okay to post it as an answer here, and if so please ping me so I can both read and up vote it. Thanks! – Gambell

@uhoh, you can follow this question. – Adjoining

@DanielWalker I've never figured that following stuff out, out but I'll look into it now, thanks! update: oh! there's a button for it right there :-) – Gambell

You should see the options "share", "edit", "follow", etc. below the tags. – Adjoining

I don't still have a good answer to this question@Gambell – Frazier

Answer 1

I deleted my previous answer - y (length n) depending on another array x, that has length k.

The current documentation suggests this is not the intention of vectorized.

Rather I think k=y.shape(2) should be interpreted as arbitrary. This means that whatever the second dimension of y is, as you pass it on to fun, fun should be able to handle that.

This means that fun should have some loop for i in y.shape(2): process column y[:, i]. Finally, return a 2D array with rhs values of shape (n, k).

The current solve_ivp documentation says

If vectorized is True, fun may be called with y of shape (n, k), where k is an integer. In this case, fun must behave such that fun(t, y)[:, i] == fun(t, y[:, i]) (i.e. each column of the returned array is the time derivative of the state corresponding with a column of y).

I think that is somewhat clearer than the original documentation (posted above).

Buskin

Interesting, for example if I have 10 non-interacting particles, and I'm solving their trajectories with six-dimensional state vectors, (x, y, z, vx, vy, vz), I can use a 2D array for y with a shape (6, 10). Here k=10. It's like running sove_ivp ten time in a row (one for each case) except that vectorize offers the possibility of sending each (non-interacting) particle to a different processor or thread or something? Is that roughly the right way to look at it? – Gambell

No, you are still solving the same initial value problem but faster if you use "Radau" or "BDF": Setting vectorized=True allows for faster finite difference approximation of the Jacobian by methods "Radau" and "BDF". – Buskin

OK got it - thank you! :-) – Gambell