Resolution systeme lorenz

Version imprimable

28/05/2015, 20h38
allmess

1 pièce(s) jointe(s)

Resolution systeme lorenz

Bnjour à tous, voilà j'ai ecris le script suivant afin de resoudre le systeme d'equations de Lorenz et 'en afficher les solutions, simplement le resultat obtenu n'est pas celui escompté, et je ne vois pas d'où proviens mon erreur..
Merci d'avance!
Pièce jointe 179096
28/05/2015, 20h53
marco056

Citation:

Envoyé par allmess

Bnjour à tous, voilà j'ai ecris le script suivant afin de resoudre le systeme d'equations de Lorenz et 'en afficher les solutions, simplement le resultat obtenu n'est pas celui escompté, et je ne vois pas d'où proviens mon erreur..
Merci d'avance!
Pièce jointe 179096

La flemme d'aller chercher le code : tu ne veux pas le mettre entre des balise de code (#) ?

Code:

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
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
 
from __future__ import absolute_import
 
import matplotlib.pyplot as pl
import numpy as np
import math
from mpl_toolkits.mplot3d import Axes3D
 
 
 
 
 
def F(x,y,z):
	return 10*(y-x)
def G(x,y,z):
	return 28*x-y-x*z
def H(x,y,z):
	return x*y-(8/3)*z
 
def Euler2(F,G,H,a,b,x0,y0,z0,n):
	x=x0
	y=y0
	z=z0
	h=(b-a)/n
	Z=[z0]
	X=[x0]
	Y=[y0]
	while x<=b:
		xi,yi,zi=x,y,z
		x=xi+h*F(xi,yi,zi)
		y=yi+h*G(xi,yi,zi)
		z=zi+h*H(xi,yi,zi)
		Z.append(z)
		X.append(x)
		Y.append(y)
		x=xi+h
	z=z0
	y=y0
	x=y0
	while x>=a:
		xi,yi,zi=x,y,z
		x=xi-h*F(xi,yi,zi)
		y=yi-h*G(xi,yi,zi)
		z=zi-h*H(xi,yi,zi)
		Z.insert(0,z)
		X.insert(0,x)
		Y.insert(0,y)
		x=xi-h
	return X,Y,Z
 
a=0
b=50
n=10000
X,Y,Z=Euler2(F,G,H,a,b,0,0.1,0,n)
#pl.plot(Y,Z)
#pl.show()
 
fig = pl.figure()
ax = fig.gca(projection='3d')
 
ax.plot(X, Z, Y, label='orbite')
ax.legend()
pl.show()

Salut.

Code:

1
2
3
4
5
6
7
8
9
10
11
def Euler2(F,G,H,a,b,x0,y0,z0,n):
    ...
    while x<=b:
        xi,yi,zi=x,y,z
        x=xi+h*F(xi,yi,zi)
        y=yi+h*G(xi,yi,zi)
        z=zi+h*H(xi,yi,zi)
        Z.append(z)
        X.append(x)
        Y.append(y)
        x=xi+h

x est modifié 2 fois dans la boucle (la seconde est de trop).
Et c'est n'est pas sur x mais sur t que porte la condition de fin. t prend n valeurs entre a et b. Une boucle for est mieux adaptée.

Code:

1
2
3
4
5
6
7
8
9
10
def Euler2(F, G, H, t1, t2, x0, y0, z0, n):
    res = [(x0, y0, z0)]
    x, y, z = res[0]
    dt = (t2 - t1) / n
    for _ in range(n):
        x += dt * F(x, y, z)
        y += dt * G(x, y, z)
        z += dt * H(x, y, z)
        res.append((x, y, z))
    return res