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