python scipy solve_ivp, почему изменение входных параметров меняет состояние системы?

Имеется следующий код:

from scipy import integrate as inter
from matplotlib import pyplot as py


def f(t,y):
    dy_dt=y[1]
    d2y_dt2=10*y[0]
    y[0]=2
    return [dy_dt,d2y_dt2]


sol=inter.solve_ivp(f,(0,1),[1,0],)

py.plot(sol.t,sol.y[0])

py.show()

Получается так, что y[0]=2 влияет на решение, хотя по идее не должно, так как функция f возвращает только производную. Почему так?


Ответы (1 шт):

Автор решения: CrazyElf

Если покопаться в исходниках SciPy (где-то здесь), то можно обнаружить, что массив y, который передаётся в вашу функцию f, используется потом для дальнейших вычислений. И ему не делается какая-то копия, вам передаётся именно тот массив, на котором производятся вычисления.

def rk_step(fun, t, y, f, h, A, B, C, K):
    ...
    y_new = y + h * np.dot(K[:-1].T, B)
    ^^^^^
    f_new = fun(t + h, y_new)
                       ^^^^^
    K[-1] = f_new

    return y_new, f_new
           ^^^^^

Если вы меняете какой-то элемент массива y внутри вашей функции, вы искажаете вычисления, производящиеся в inter.solve_ivp, что и наблюдается наглядно на графике.

→ Ссылка