Как установить фиксированный размер шага с помощью scipy.integrate?

Таблицу Мясника для метода Дорманда-Принца RK45 довольно легко запрограммировать.

0
1/5   |  1/5
3/10  |  3/40        9/40
4/5   |  44/45       −56/15        32/9
8/9   |  19372/6561  −25360/2187   64448/6561   −212/729
1     |  9017/3168   −355/33       46732/5247   49/176     −5103/18656
1     |  35/384           0        500/1113     125/192    −2187/6784     11/84     
-----------------------------------------------------------------------------------------
      |  35/384           0        500/1113     125/192    −2187/6784     11/84     0
      |  5179/57600       0        7571/16695   393/640    −92097/339200  187/2100  1/40

сначала в функции для одношагового импорта numpy как np

def DoPri45Step(f,t,x,h):

    k1 = f(t,x)
    k2 = f(t + 1./5*h, x + h*(1./5*k1) )
    k3 = f(t + 3./10*h, x + h*(3./40*k1 + 9./40*k2) )
    k4 = f(t + 4./5*h, x + h*(44./45*k1 - 56./15*k2 + 32./9*k3) )
    k5 = f(t + 8./9*h, x + h*(19372./6561*k1 - 25360./2187*k2 + 64448./6561*k3 - 212./729*k4) )
    k6 = f(t + h, x + h*(9017./3168*k1 - 355./33*k2 + 46732./5247*k3 + 49./176*k4 - 5103./18656*k5) )

    v5 = 35./384*k1 + 500./1113*k3 + 125./192*k4 - 2187./6784*k5 + 11./84*k6
    k7 = f(t + h, x + h*v5)
    v4 = 5179./57600*k1 + 7571./16695*k3 + 393./640*k4 - 92097./339200*k5 + 187./2100*k6 + 1./40*k7;

    return v4,v5

а затем в стандартной петле с фиксированным шагом

def DoPri45integrate(f, t, x0):
    N = len(t)
    x = [x0]
    for k in range(N-1):
        v4, v5 = DoPri45Step(f,t[k],x[k],t[k+1]-t[k])
        x.append(x[k] + (t[k+1]-t[k])*v5)
    return np.array(x)

Затем протестируйте его на каком-нибудь игрушечном примере с известным точным решением y(t)=sin(t)

def mms_ode(t,y): return np.array([ y[1], sin(sin(t))-sin(t)-sin(y[0]) ])
mms_x0 = [0.0, 1.0]

и построить график ошибки в масштабе h^5

for h in [0.2, 0.1, 0.08, 0.05, 0.01][::-1]:
    t = np.arange(0,20,h);
    y = DoPri45integrate(mms_ode,t,mms_x0)
    plt.plot(t, (y[:,0]-np.sin(t))/h**5, 'o', ms=3, label = "h=%.4f"%h);
plt.grid(); plt.legend(); plt.show()  

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

Посмотрев на реализацию шага, вы обнаружите, что лучшее, что вы можете сделать, - это контролировать начальный размер шага (в пределах, установленных минимальным и максимальным размером шага), установив атрибут h_abs перед вызовом RK45.step:

In [27]: rk = RK45(lambda t, y: t, 0, [0], 1e6)

In [28]: rk.h_abs = 30

In [29]: rk.step()

In [30]: rk.step_size
Out[30]: 30.0

Scipy.integrate обычно используется с методом изменяемого шага, контролируя TOL (ошибку одного шага) при численном интегрировании. TOL обычно вычисляется путем проверки другим численным методом. Например, RK45 использует метод Рунге-Кутты 5-го порядка для проверки TOL метода Рунге-Кутта 4-го порядка для определения шага интегрирования.

Следовательно, если вам необходимо интегрировать ODE с фиксированным шагом, просто отключите проверку TOL, установив atol, rtol с довольно большой константой. Например, как форма:

solve_ivp(your function, t_span=[0, 10], y0=..., method="RK45", max_step=0.01, atol = 1, rtol = 1)

Проверка TOL установлена ​​настолько большой, что шагом интегрирования будет выбранный вами max_step.

Вы сказали, что хотите поведение с фиксированным временным шагом, а не просто с фиксированным временным шагом эволюции. Следовательно, вам придется проделать свой путь через это, если вы не хотите заново реализовывать решатель самостоятельно. Просто установите допуски интеграции atol и rtol на 1e90, а max_step и first_step на значение dt временного шага, который вы хотите использовать. Таким образом, расчетная ошибка интегрирования всегда будет очень маленькой, что заставит решатель не сокращать временной шаг динамически.

Однако используйте этот прием только с ЯВНЫМИ алгоритмами (RK23, RK45, DOP853)! Неявные алгоритмы от resolve_ivp (Radau, BDF, возможно, LSODA также) регулируют точность нелинейного решателя Ньютона в соответствии с atol и rtol, поэтому вы можете получить решение, которое не имеет никакого смысла ...

Если вас интересует размер шага исправления с учетом данных, я настоятельно рекомендую вам использовать функцию scipy.integrate.solve_ivp и ее аргумент t_eval.

Эта функция объединяет все scipy.integrate решатели од в одну функцию, поэтому вам нужно выбрать метод, присвоив значение его method аргументу. К счастью, метод по умолчанию - RK45, так что вам не нужно беспокоиться об этом.

Что вам более интересно, так это аргумент t_eval, в котором вы должны указать плоский массив. Функция производит выборку кривой решения при t_eval значениях и возвращает только эти точки. Поэтому, если вам нужна единообразная выборка по размеру шага, просто укажите аргументу t_eval следующее: numpy.linspace(t0, tf, samplingResolution), где t0 - начало, а tf - конец моделирования. Таким образом, вы можете иметь единообразную выборку без необходимости прибегать к фиксированному размеру шага, который вызывает нестабильность для некоторых ODE.

Я предлагаю написать свою собственную программу с фиксированным шагом rk4 на py. Есть много примеров из Интернета, которые могут вам помочь. Это гарантирует, что вы точно знаете, как вычисляется каждое значение. Кроме того, обычно не будет вычислений 0/0, и если да, то их будет легко отследить и побудить еще раз взглянуть на решаемую оду.