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

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

0 голосов
/ 02 февраля 2019

Я ищу способ установить фиксированный размер шага для решения моей проблемы начальных значений методом Рунге-Кутты в Python.Соответственно, как я могу сказать scipy.integrate.RK45 постоянно обновлять (размер шага) процедуры его интеграции?

Большое спасибо.

Ответы [ 2 ]

0 голосов
/ 03 февраля 2019

Довольно просто закодировать таблицу Мясника для метода Dormand-Prince 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

сначала в функции за один шаг 

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 = np.asarray(N*[x0]);
    for k in range(N-1):
        v4, v5 = DoPri45Step(f,t[k],x[k],t[k+1]-t[k])
        x[k+1] = x[k] + (t[k+1]-t[k])*v5
    return 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, поскольку графики коэффициентов ошибок сближаются.

enter image description here

0 голосов
/ 02 февраля 2019

Посмотрев на реализацию шага , вы обнаружите, что лучшее, что вы можете сделать, это контролировать начальный размер шага (в пределах, установленных минимальным и максимальным размером шага), установиватрибут 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
...