Comment définir une taille de pas fixe avec scipy.integrate?

En consultant l'implémentation de l'étape , vous constaterez que le mieux que vous puissiez faire est de contrôler la taille initiale du pas (dans les limites fixées par la taille minimale et maximale du pas) en définissant l'attribut h_abs avant d'appeler 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

Il est assez facile de coder le tableau Butcher pour la méthode Dormand-Prince RK45.

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()  

d'abord dans une fonction pour une seule étape import numpy as np

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

puis dans une boucle standard à pas fixe

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)

Puis testez-le pour un exemple de jouet avec la solution exacte connue y(t)=sin(t)

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

et tracez l'erreur mise à l'échelle par h ^ 5 code >

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

pour obtenir la confirmation qu'il s'agit bien d'une méthode d'ordre 5, car les graphiques des coefficients d'erreur se rapprochent.

Scipy.integrate est généralement utilisé avec une méthode d'étape modifiable en contrôlant le TOL (erreur d'une étape) tout en intégrant numériquement. Le TOL est généralement calculé en vérifiant avec une autre méthode numérique. Par exemple, RK45 utilise le 5ème ordre Runge-Kutta pour vérifier le TOL de la 4ème méthode Runge-Kutta d'ordre pour déterminer l'étape d'intégration.

Par conséquent, si vous devez intégrer des ODE avec un pas fixe, désactivez simplement le contrôle TOL en définissant atol, rtol avec une constante assez grande. Par exemple, comme le formulaire:

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

La vérification TOL est définie pour être si grande que l'étape d'intégration correspond au max_step que vous choisissez.

Si vous êtes intéressé par la taille de l'étape de correction des données, je vous recommande vivement d'utiliser la fonction scipy.integrate.solve_ivp et son argument t_eval .

Cette fonction regroupe tous les solveurs ode scipy.integrate en une seule fonction, vous devez donc choisir la méthode en donnant une valeur à son argument method . Heureusement, la méthode par défaut est le RK45, vous n'avez donc pas à vous en préoccuper.

Ce qui est plus intéressant pour vous, c'est l'argument t_eval , où vous devez donner un tableau plat. La fonction échantillonne la courbe de solution aux valeurs de t_eval et ne renvoie que ces points. Donc, si vous voulez un échantillonnage uniforme par la taille du pas, donnez simplement à l'argument t_eval ce qui suit: numpy.linspace (t0, tf, samplingResolution) , où t0 est le début et tf est la fin de la simulation. Ainsi, vous pouvez avoir un échantillonnage uniforme sans avoir à recourir à une taille de pas fixe qui provoque une instabilité pour certains ODE.

Vous avez dit que vous souhaitiez un comportement de pas de temps fixe, pas seulement un pas de temps d'évolution fixe. Par conséquent, vous devez "pirater" votre chemin si vous ne voulez pas réimplémenter le solveur vous-même. Il vous suffit de définir les tolérances d'intégration atol et rtol sur 1e90, et max_step et first_step sur la valeur dt du pas de temps que vous souhaitez utiliser . De cette façon, l'erreur d'intégration estimée sera toujours très faible, ce qui incitera le solveur à ne pas réduire le pas de temps de manière dynamique.

Cependant, n'utilisez cette astuce qu'avec les algortithmes EXPLICIT (RK23, RK45, DOP853)! Les algorithmes implicites de "résoudre_ivp" (Radau, BDF, peut-être aussi LSODA) ajustent la précision du solveur de Newton non linéaire en fonction de atol et rtol, donc vous pourriez finir par avoir une solution qui n'a aucun sens ...

Je suggère d'écrire votre propre programme à étapes fixes rk4 dans py. Il existe de nombreux exemples Internet pour vous aider. Cela garantit que vous savez précisément comment chaque valeur est calculée. De plus, il n'y aura normalement pas de calculs 0/0 et si c'est le cas, ils seront faciles à tracer et inciteront à un autre regard sur l'ode en cours de résolution.