Systèmes différentiels : Méthode de Runge Kutta ordre 4 - équation de Lorenz / Navier-Stockes & attracteur de Lorenz

Chaka · le 25 avril 2021

Bonjour,
Dans le cadre d'un projet d'étude (3ème année de Licence), je voudrais mettre en place, sur Python, la méthode de RK4 pour résoudre les équations de Lorentz numériquement.
J'ai réussi à le faire sans trop de difficulté avec les méthodes de Newton, mais pour RK4, je bloque !! Le calcul des intermédiaires m'est compliqué.
Je vous mets ci-joint mes scripts Python.
Merci par avance pour vos retours.

Chaka

PS: ce post à déjà été émis dans le forum Mathématiques.

Projet P.py

julesx · le 25 avril 2021

Bonsoir,

Ton post en maths a été verrouillé, mais c'est accessoire. Par contre, à mon avis, pour que tu aies une aide, je pense qu'il faudrait que tu réduises ton script à la seule partie qui pose problème. Si c'est possible évidemment, je ne suis pas assez compétent dans ce domaine pour en juger.

Chaka · le 25 avril 2021

Je vais reposter mon script de la partie qui m'intéresse dans ce cas . Merci

Chaka · le 26 avril 2021

Lorenz avec la méthode d'Euler :

def lorenz(x, y, z, s=10, r=28, b=2.667):

"""

Entrée:

x, y, z: point dans R3 auquel nous nous interessons

s, r, b: paramètres de l'attracteur de Lorenz

Sortie:

dérivées partiels aux points x, y, z liées aux équations de Lorenz

"""

x_dot = s*(y - x)

y_dot = r*x - y - x*z

z_dot = x*y - b*z

return x_dot, y_dot, z_dot

pas = 0.01

num_iter = 10000

# +1 pour prendre les CI

x_sol = np.empty(num_iter + 1)

y_sol = np.empty(num_iter + 1)

z_sol = np.empty(num_iter + 1)

# CI : Conditions initiales

x_sol[0], y_sol[0], z_sol[0] = (0., 1., 1.05)

# Euler

for i in range(num_iter):

x_dot, y_dot, z_dot = lorenz(x_sol[i], y_sol[i], z_sol[i])

x_sol[i + 1] = x_sol[i] + (x_dot * pas)

y_sol[i + 1] = y_sol[i] + (y_dot * pas)

z_sol[i + 1] = z_sol[i] + (z_dot * pas)

# Plot

ax = plt.figure().add_subplot(projection='3d')

ax.plot(x_sol, y_sol, z_sol, lw=0.5)

ax.set_xlabel("X")

ax.set_ylabel("Y")

ax.set_zlabel("Z")

ax.set_title("Attracteur de Lorenz")

plt.show()

Runge Kutta d'odre 4 :

def rk4(f,t0,tf,y0,N):

y=[y0]

t=np.linspace(t0,tf,N+1)

pas=(tf-t0)/N

#Initialisation

k1=[f(y[i],t[i])]

k2=[y0+(pas/2)*f(t0,y0)]

k3=[y0+(pas/2)*f(t0+pas/2,k1[0])]

k4=[y0+pas*f(t0+(pas/2),k2[0])]

for i in range(N):

y.append(y[i]+(pas/6)*(f(t[i],y[i])+2*f(t[i]+pas/2,k2[i])+2*f(t[i]+(pas/2),k3[i-1])+f(t[i+1],k4[i])))

k2.append(y[i+1]+(pas/2)*f(t[i+1],y[i+1]))

k3.append(y[i+1]+(pas/2)*f(t[i+1]+(pas/2),k2[i+1]))

return t,y

julesx · le 26 avril 2021

Re-bonjour,

Je ne sais pas si ça va pouvoir te servir. J'ai trouvé ça sur la toile à l'adresse

https://fr.mathworks.com/matlabcentral/answers/360751-solving-lorenz-attractor-equations-using-runge-kutta-rk4-method

 clc;
 clear all;
t(1)=0;  %initializing x,y,z,t
x(1)=1;
y(1)=1;
z(1)=1;
sigma=10;   %value of constants
rho=28;
beta=(8.0/3.0);
h=0.01;   %step size
t=0:h:20;
f=@(t,x,y,z) sigma*(y-x);  %ode 
g=@(t,x,y,z) x*rho-x.*z-y;
p=@(t,x,y,z) x.*y-beta*z;
for i=1:(length(t)-1) %loop
    k1=f(t(i),x(i),y(i),z(i));
    l1=g(t(i),x(i),y(i),z(i));
    m1=p(t(i),x(i),y(i),z(i));
      k2=f(t(i)+h/2,(x(i)+0.5*k1*h),(y(i)+(0.5*l1*h)),(z(i)+(0.5*m1*h)));     
      l2=g(t(i)+h/2,(x(i)+0.5*k1*h),(y(i)+(0.5*l1*h)),(z(i)+(0.5*m1*h)));
      m2=p(t(i)+h/2,(x(i)+0.5*k1*h),(y(i)+(0.5*l1*h)),(z(i)+(0.5*m1*h)));
      k3=f(t(i)+h/2,(x(i)+0.5*k2*h),(y(i)+(0.5*l2*h)),(z(i)+(0.5*m2*h)));
      l3=g(t(i)+h/2,(x(i)+0.5*k2*h),(y(i)+(0.5*l2*h)),(z(i)+(0.5*m2*h)));
      m3=p(t(i)+h/2,(x(i)+0.5*k2*h),(y(i)+(0.5*l2*h)),(z(i)+(0.5*m2*h)));
      k4=f(t(i)+h,(x(i)+k3*h),(y(i)+l3*h),(z(i)+m3*h));
      l4=g(t(i)+h,(x(i)+k3*h),(y(i)+l3*h),(z(i)+m3*h));
      m4=p(t(i)+h,(x(i)+k3*h),(y(i)+l3*h),(z(i)+m3*h));
      x(i+1) = x(i) + h*(k1 +2*k2  +2*k3   +k4)/6; %final equations
      y(i+1) = y(i) + h*(l1  +2*l2   +2*l3    +l4)/6;
      z(i+1) = z(i) + h*(m1+2*m2 +2*m3  +m4)/6;
end
plot3(x,y,z)

C'est du Matlab, ça devrait pouvoir se transcrire. Il n'y a que les lignes

f=@(t,x,y,z) sigma*(y-x); %ode

g=@(t,x,y,z) x*rho-x.*z-y;

p=@(t,x,y,z) x.*y-beta*z;

qui sont à "décoder", mais ça correspond peut-être à la notion de def de Python.

Mais ne m'en demande pas plus, je n'ai jamais utilisé Matlab autrement qu'avec la sous-couche Simulink.

Chaka · le 26 avril 2021

Merci beaucoup julesx !!
Je viens de lire le dernier post de ce forum et cela à l'air de répondre à ma question.
Entre temps j'avais un peu avancé et étais en bonne voie d'après ce que je viens de lire.
Je t'en remercie !

Je posterai le code final pour ceux que cela intéresse
Bonne soirée

Chaka · le 26 avril 2021

Je viens de finir mon script cependant je rencontre un problème par rapport aux résultats que j'obtiens : cf pièces jointes et scripts.
La première est la figure que j'obtiens avec RK4 (convergence vers l'origine), la seconde avec la méthode de Newton (aspect chaotique). (avec les mêmes constantes et les mêmes conditions initiales)
Trouvez-vous une explication à ces différences ? Si non, problèmes dans les scripts ? :

## Euler NS-Lorenz

def lorenz(x, y, z, s=10, r=28, b=2.667):

"""

Entrée:

x, y, z: point dans R3 auquel nous nous interessons

s, r, b: paramètres de l'attracteur de Lorenz

Sortie:

dérivées partiels aux points x, y, z liées aux équations de Lorenz

"""

x_dot = s*(y - x)

y_dot = r*x - y - x*z

z_dot = x*y - b*z

return x_dot, y_dot, z_dot

pas = 0.01

num_iter = 10000

# +1 pour prendre les CI

x_sol = np.empty(num_iter + 1)

y_sol = np.empty(num_iter + 1)

z_sol = np.empty(num_iter + 1)

# CI : Conditions initiales

x_sol[0], y_sol[0], z_sol[0] = (0., 1., 1.05)

# Euler

for i in range(num_iter):

x_dot, y_dot, z_dot = lorenz(x_sol[i], y_sol[i], z_sol[i])

x_sol[i + 1] = x_sol[i] + (x_dot * pas)

y_sol[i + 1] = y_sol[i] + (y_dot * pas)

z_sol[i + 1] = z_sol[i] + (z_dot * pas)

# Plot

ax = plt.figure().add_subplot(projection='3d')

ax.plot(x_sol, y_sol, z_sol, lw=0.5)

ax.set_xlabel("X")

ax.set_ylabel("Y")

ax.set_zlabel("Z")

ax.set_title("Attracteur de Lorenz")

plt.show()

## RK4 NS-Lorenz

s,r,b=10,28,2.667

#Lorenz

def f(t,x,y,z):

return s*(y-x)

def g(t,x,y,z):

return r*x-y-x*z

def h(t,x,y,z):

return x*y-b*z

# Initialisation

x,y,z=[0],[1],[1.05]

pas=0.01

num_iter=10000 #t=pas*num_iter

t=np.linspace(0,pas*num_iter,num_iter+1)

# RK4

for i in range (num_iter):

k1=pas*f(t[i],x[i],y[i],z[i])

l1=pas*g(t[i],x[i],y[i],z[i])

m1=pas*f(t[i],x[i],y[i],z[i])

k2=f(t[i]+pas/2,(x[i]+k1*pas/2),(y[i]+(l1*pas/2)),(z[i]+(m1*pas/2)))

l2=g(t[i]+pas/2,(x[i]+k1*pas/2),(y[i]+(l1*pas/2)),(z[i]+(m1*pas/2)))

m2=h(t[i]+pas/2,(x[i]+k1*pas/2),(y[i]+(l1*pas/2)),(z[i]+(m1*pas/2)))

k3=f(t[i]+pas/2,(x[i]+k1*pas/2),(y[i]+(l1*pas/2)),(z[i]+(m1*pas/2)))

l3=g(t[i]+pas/2,(x[i]+k1*pas/2),(y[i]+(l1*pas/2)),(z[i]+(m1*pas/2)))

m3=h(t[i]+pas/2,(x[i]+k1*pas/2),(y[i]+(l1*pas/2)),(z[i]+(m1*pas/2)))

k4=f(t[i]+pas,(x[i]+k3*pas),(y[i]+l3*pas),(z[i]+m3*pas))

l4=g(t[i]+pas,(x[i]+k3*pas),(y[i]+l3*pas),(z[i]+m3*pas))

m4=h(t[i]+pas,(x[i]+k3*pas),(y[i]+l3*pas),(z[i]+m3*pas))

x.append(x[i]+pas*(k1+2*k2+2*k3+k4)/6)

y.append(y[i]+pas*(l1+2*l2+2*l3+l4)/6)

z.append(z[i]+pas*(k1+2*k2+2*k3+k4)/6)

# Plot

ax = plt.figure().add_subplot(projection='3d')

ax.plot(x, y, z, lw=0.5)

ax.set_xlabel("X")

ax.set_ylabel("Y")

ax.set_zlabel("Z")

ax.set_title("Attracteur de Lorenz")

plt.show()

Modifié le 26 avril 2021 par Chaka

Chaka · le 27 avril 2021

Il y a 10 heures, Chaka a dit :

Je viens de finir mon script cependant je rencontre un problème par rapport aux résultats que j'obtiens : cf pièces jointes et scripts.
La première est la figure que j'obtiens avec RK4 (convergence vers l'origine), la seconde avec la méthode de Newton (aspect chaotique). (avec les mêmes constantes et les mêmes conditions initiales)
Trouvez-vous une explication à ces différences ? Si non, problèmes dans les scripts ? :

## Euler NS-Lorenz

def lorenz(x, y, z, s=10, r=28, b=2.667):

"""

Entrée:

x, y, z: point dans R3 auquel nous nous interessons

s, r, b: paramètres de l'attracteur de Lorenz

Sortie:

dérivées partiels aux points x, y, z liées aux équations de Lorenz

"""

x_dot = s*(y - x)

y_dot = r*x - y - x*z

z_dot = x*y - b*z

return x_dot, y_dot, z_dot

pas = 0.01

num_iter = 10000

# +1 pour prendre les CI

x_sol = np.empty(num_iter + 1)

y_sol = np.empty(num_iter + 1)

z_sol = np.empty(num_iter + 1)

# CI : Conditions initiales

x_sol[0], y_sol[0], z_sol[0] = (0., 1., 1.05)

# Euler

for i in range(num_iter):

x_dot, y_dot, z_dot = lorenz(x_sol[i], y_sol[i], z_sol[i])

x_sol[i + 1] = x_sol[i] + (x_dot * pas)

y_sol[i + 1] = y_sol[i] + (y_dot * pas)

z_sol[i + 1] = z_sol[i] + (z_dot * pas)

# Plot

ax = plt.figure().add_subplot(projection='3d')

ax.plot(x_sol, y_sol, z_sol, lw=0.5)

ax.set_xlabel("X")

ax.set_ylabel("Y")

ax.set_zlabel("Z")

ax.set_title("Attracteur de Lorenz")

plt.show()

## RK4 NS-Lorenz

s,r,b=10,28,2.667

#Lorenz

def f(t,x,y,z):

return s*(y-x)

def g(t,x,y,z):

return r*x-y-x*z

def h(t,x,y,z):

return x*y-b*z

# Initialisation

x,y,z=[0],[1],[1.05]

pas=0.01

num_iter=10000 #t=pas*num_iter

t=np.linspace(0,pas*num_iter,num_iter+1)

# RK4

for i in range (num_iter):

k1=pas*f(t[i],x[i],y[i],z[i])

l1=pas*g(t[i],x[i],y[i],z[i])

m1=pas*f(t[i],x[i],y[i],z[i])

k2=f(t[i]+pas/2,(x[i]+k1*pas/2),(y[i]+(l1*pas/2)),(z[i]+(m1*pas/2)))

l2=g(t[i]+pas/2,(x[i]+k1*pas/2),(y[i]+(l1*pas/2)),(z[i]+(m1*pas/2)))

m2=h(t[i]+pas/2,(x[i]+k1*pas/2),(y[i]+(l1*pas/2)),(z[i]+(m1*pas/2)))

k3=f(t[i]+pas/2,(x[i]+k1*pas/2),(y[i]+(l1*pas/2)),(z[i]+(m1*pas/2)))

l3=g(t[i]+pas/2,(x[i]+k1*pas/2),(y[i]+(l1*pas/2)),(z[i]+(m1*pas/2)))

m3=h(t[i]+pas/2,(x[i]+k1*pas/2),(y[i]+(l1*pas/2)),(z[i]+(m1*pas/2)))

k4=f(t[i]+pas,(x[i]+k3*pas),(y[i]+l3*pas),(z[i]+m3*pas))

l4=g(t[i]+pas,(x[i]+k3*pas),(y[i]+l3*pas),(z[i]+m3*pas))

m4=h(t[i]+pas,(x[i]+k3*pas),(y[i]+l3*pas),(z[i]+m3*pas))

x.append(x[i]+pas*(k1+2*k2+2*k3+k4)/6)

y.append(y[i]+pas*(l1+2*l2+2*l3+l4)/6)

z.append(z[i]+pas*(k1+2*k2+2*k3+k4)/6)

# Plot

ax = plt.figure().add_subplot(projection='3d')

ax.plot(x, y, z, lw=0.5)

ax.set_xlabel("X")

ax.set_ylabel("Y")

ax.set_zlabel("Z")

ax.set_title("Attracteur de Lorenz")

plt.show()

Bonjour,
Je viens de trouver mon erreur! dans la ligne z.append(...), il faut remplacer les k1,k2,... par m1,m2,... (j'ai voulu faire un copié collé pour gagner du temps mais finalement cela m'en a fait perdre)
Bonne journée !

Chaka

julesx · le 27 avril 2021

Bonjour,

Effectivement, avec la rectification, on obtient bien les mêmes résultats.

Juste une remarque, à l'avenir, pour des scripts, il vaut mieux les publier en utilisant la balise code < > car ainsi, les indentations sont parfaitement respectées.

Bonne journée également.

Daddy06 · le 11 septembre 2021

Le 26/04/2021 à 10:08, Chaka a dit :

Lorenz avec la méthode d'Euler :

def lorenz(x, y, z, s=10, r=28, b=2.667):

"""

Entrée:

x, y, z: point dans R3 auquel nous nous interessons

s, r, b: paramètres de l'attracteur de Lorenz

Sortie:

dérivées partiels aux points x, y, z liées aux équations de Lorenz

"""

x_dot = s*(y - x)

y_dot = r*x - y - x*z

z_dot = x*y - b*z

return x_dot, y_dot, z_dot

pas = 0.01

num_iter = 10000

# +1 pour prendre les CI

x_sol = np.empty(num_iter + 1)

y_sol = np.empty(num_iter + 1)

z_sol = np.empty(num_iter + 1)

# CI : Conditions initiales

x_sol[0], y_sol[0], z_sol[0] = (0., 1., 1.05)

# Euler

for i in range(num_iter):

x_dot, y_dot, z_dot = lorenz(x_sol[i], y_sol[i], z_sol[i])

x_sol[i + 1] = x_sol[i] + (x_dot * pas)

y_sol[i + 1] = y_sol[i] + (y_dot * pas)

z_sol[i + 1] = z_sol[i] + (z_dot * pas)

# Plot

ax = plt.figure().add_subplot(projection='3d')

ax.plot(x_sol, y_sol, z_sol, lw=0.5)

ax.set_xlabel("X")

ax.set_ylabel("Y")

ax.set_zlabel("Z")

ax.set_title("Attracteur de Lorenz")

plt.show()

Runge Kutta d'odre 4 :
[url=https://www.rachat-de-credit-simulation.com/rachat-de-credit]Rachat de crédit et trésorerie[/url]

def rk4(f,t0,tf,y0,N):

y=[y0]

t=np.linspace(t0,tf,N+1)

pas=(tf-t0)/N

#Initialisation

k1=[f(y[i],t[i])]

k2=[y0+(pas/2)*f(t0,y0)]

k3=[y0+(pas/2)*f(t0+pas/2,k1[0])]

k4=[y0+pas*f(t0+(pas/2),k2[0])]

for i in range(N):

y.append(y[i]+(pas/6)*(f(t[i],y[i])+2*f(t[i]+pas/2,k2[i])+2*f(t[i]+(pas/2),k3[i-1])+f(t[i+1],k4[i])))

k2.append(y[i+1]+(pas/2)*f(t[i+1],y[i+1]))

k3.append(y[i+1]+(pas/2)*f(t[i+1]+(pas/2),k2[i+1]))

return t,y

Oh lala! C'est compliqué, j'en ai mal à la tête 😕

julesx · le 11 septembre 2021

Bonjour,

Mais tu as regardé ce post car tu cherchais une information ou c'était uniquement pour te faire peur ?

Connexion

Systèmes différentiels : Méthode de Runge Kutta ordre 4 - équation de Lorenz / Navier-Stockes & attracteur de Lorenz

Messages recommandés

Chaka

Lien vers le commentaire

Partager sur d’autres sites

julesx

Lien vers le commentaire

Partager sur d’autres sites

Chaka

Lien vers le commentaire

Partager sur d’autres sites

Chaka

Lien vers le commentaire

Partager sur d’autres sites

julesx

Lien vers le commentaire

Partager sur d’autres sites

Chaka

Lien vers le commentaire

Partager sur d’autres sites

Chaka

Lien vers le commentaire

Partager sur d’autres sites

Chaka

Lien vers le commentaire

Partager sur d’autres sites

julesx

Lien vers le commentaire

Partager sur d’autres sites

Daddy06

Lien vers le commentaire

Partager sur d’autres sites

julesx

Lien vers le commentaire

Partager sur d’autres sites

Rejoindre la conversation

Naviguer

Activité