Méthode de quatrième ordre de Runge-Kutta

Question

J'applique la méthode de quatrième ordre de Runge-Kutta pour le système de deux équations.

h est le nombre de segments, de sorte que T/h est pas.

def cauchy(f1, f2, x10, x20, T, h): 
    x1 = [x10] 
    x2 = [x20] 

    for i in range(1, h): 
     k11 = f1((i-1)*T/h, x1[-1], x2[-1]) 
     k12 = f2((i-1)*T/h, x1[-1], x2[-1]) 
     k21 = f1((i-1)*T/h + T/h/2, x1[-1] + T/h/2*k11, x2[-1] + T/h/2*k12) 
     k22 = f2((i-1)*T/h + T/h/2, x1[-1] + T/h/2*k11, x2[-1] + T/h/2*k12) 
     k31 = f1((i-1)*T/h + T/h/2, x1[-1] + T/h/2*k21, x2[-1] + T/h/2*k22) 
     k32 = f2((i-1)*T/h + T/h/2, x1[-1] + T/h/2*k21, x2[-1] + T/h/2*k22) 
     k41 = f1((i-1)*T/h + T/h, x1[-1] + T/h*k31, x2[-1] + T/h*k32) 
     k42 = f2((i-1)*T/h + T/h, x1[-1] + T/h*k31, x2[-1] + T/h*k32) 

     x1.append(x1[-1] + T/h/6*(k11 + 2*k21 + 2*k31 + k41)) 
     x2.append(x2[-1] + T/h/6*(k12 + 2*k22 + 2*k32 + k42)) 

    return x1, x2 

Ensuite, je teste sur ce système:

def f1(t, x1, x2): 
    return x2 

def f2(t, x1, x2): 
    return -x1 

def true_x1(t): 
    return np.cos(t) + np.sin(t) 

def true_x2(t): 
    return np.cos(t) - np.sin(t) 

Il semble fonctionner très bien (je l'ai testé aussi avec des valeurs initiales différentes et fonctions différentes: tous les travaux fin):

x10 = 1 
x20 = 1 
T = 1 
h = 10 

x1, x2 = cauchy(f1, f2, x10, x20, T, h) 
t = np.linspace(0, T, h) 

plt.xlabel('t') 
plt.ylabel('x1') 
plt.plot(t, true_x1(t), "blue", label="true_x1") 
plt.plot(t, x1, "red", label="approximation_x1") 
plt.legend(bbox_to_anchor=(0.97, 0.27)) 
plt.show() 

plt.xlabel('t') 
plt.ylabel('x2') 
plt.plot(t, true_x2(t), "blue", label="true_x2") 
plt.plot(t, x2, "red", label="approximation_x2") 
plt.legend(bbox_to_anchor=(0.97, 0.97)) 
plt.show() 

Je veux vérifier si l'erreur est de l'ordre de O(step^4), donc je réduis pas et erreur de calcul comme ceci:

step = [] 
x1_error = [] 
x2_error = [] 
for segm in reversed(range(10, 1000)): 
    x1, x2 = cauchy(f1, f2, x10, x20, T, segm) 
    t = np.linspace(0, T, segm) 
    step.append(1/segm) 
    x1_error.append(np.linalg.norm(x1 - true_x1(t), np.inf)) 
    x2_error.append(np.linalg.norm(x2 - true_x2(t), np.inf)) 

Et je reçois ceci:

plt.plot(step, x1_error, label="x1_error") 
plt.plot(step, x2_error, label="x2_error") 
plt.legend() 

Donc, l'erreur est linéaire à partir de l'étape. C'est vraiment étrange, car il est censé être de l'ordre de O(step^4). Quelqu'un peut-il me dire ce que j'ai fait de mal?

Answer 1

for i in range(1, h): 

Cela itérer 1-h-1. Comme la dernière étape est manquante, la différence de x[h-1] à l'instant T-T/h à la solution exacte à l'instant T est O(T/h).

Ainsi, l'utilisation

for i in range(1,h+1): 

pour h pas de i-1-i ou

for i in range(h): 

pour h pas i-i+1.

---

De plus, np.linspace(0,1,4) va produire 4 nombres équidistants où le premier est 0 et le dernier est 1, ce qui

array([ 0.  , 0.33333333, 0.66666667, 1.  ]) 

qui est sans doute pas ce que vous attendiez. Ainsi, avec la correction ci-dessus, utiliser

t = np.linspace(0, T, segm+1) 

pour utiliser les mêmes points de temps dans les deux calculs.

---

Il serait plus facile de suivre votre code si vous utilisez les lettres dans leur sens habituel, où h ou dt est la taille de l'étape et N est le nombre d'étapes. Puis définissez avant la boucle h=T/N ou dt=T/N pour éviter l'utilisation répétée de T/N dans les appels de fonction.