Discussion :

[INoctal]

Bonjour,

Dans le cadre d'un TP, je dois programmer une procédure permettant de résoudre des équations différentielles avec la méthode RK4.
J'y suis bien arrivé pour une équation du premier ordre mais je bloque quand il s'agit de traiter les équations du second ..

Voici le code que j'ai écrit, (pour l'équation d²x/dt² + x = 0):
(J'ai posé v=dx/dt)

import numpy as np
import matplotlib.pyplot as plt

ti=0
tf=5
h=0.1

t=list(np.arange(ti,tf,h))

n=len(t)

def F1(x,v):
return(v)

def F2(x,v):
return -x

def RK4(F1,F2):
x=1
v=2
X=[x]

for i in range(n-1):
xk1=F1(x,v)
vk1=F2(x,v)

xk2=F1(x,x+h/2*xk1)
vk2=F2(x,v+h/2*vk1)

xk3=F1(x,x+h/2*xk2)
vk3=F2(x,v+h/2*vk2)

xk4=F1(x,x+h*xk3)
vk4=F2(x,v+h*vk3)

x=x+h/6*(xk1+2*xk2+2*xk3+xk4)
v=v+h/6*(vk1+2*vk2+2*vk3+vk4)

X.append(x)

print(X)
plt.plot(t,X)
plt.show()

Le code ne marche pas et voila ce que j'obtient:




Je vous remercie par avance,

INoctal

[INoctal]

J'ai finalement trouvé la solution en réécrivant totalement la procédure.

Pour ceux que ça pourrait intéresser, voila mon code final:

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
h=0.1
t=list(np.arange(0,100,h))
n=len(t)
 
U=np.array([w,dw])
 
def f1(w,dw):
    return dw
def f2(w,dw):  # l'équation ici est d²w/dt + w = 0
    return -w
 
def f(U):
    return np.array([f1,f2]) 
 
def RK4(f):
    w=0
    dw=1
    W=[w]
    for k in range(n-1):
        k11=f1(w,dw)
        k12=f2(w,dw)
 
        k21=f1(w+h/2*k11,dw+h/2*k12)
        k22=f2(w+h/2*k11,dw+h/2*k12)
 
        k31=f1(w+h/2*k21,dw+h/2*k22)
        k32=f2(w+h/2*k21,dw+h/2*k22)
 
        k41=f1(w+h*k31, dw+h*k32)
        k42=f2(w+h*k31, dw+h*k32)
 
        w=w+h*(k11+2*k21+2*k31+k41)/6
        dw=dw+h*(k12+2*k22+2*k32+k42)/6
        W.append(w)
 
    plt.plot(t,W)
    plt.show