# Demo entry 6357192

lorenz

Submitted by anonymous on Apr 21, 2017 at 10:59
Language: Python 3. Code size: 1.2 kB.

```# -*- coding: utf-8 -*-
"""
Created on Fri Apr 21 16:58:51 2017

@author: yangd
"""

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint
from mpl_toolkits.mplot3d import Axes3D
import math

#赋值区
T_end=1e4
T_step=1e-3
N_delta=1/8
delta=1/N_delta

# 给出位置矢量w，和三个参数sigma, r, b计算出dx/dt, dy/dt, dz/dt的值
def lorenz(w, t, sigma, r, b):
x, y, z = w
return np.array([sigma*(y-x), x*(r-z)-y, x*y-b*z])

# 创建时间点
t = np.arange(0, T_end, T_step)

#初始条件
track = odeint(lorenz, (0.01,0.01,0.01), t, args=(10.0,28,8.0/3.0))

# 3D作图
fig = plt.figure()
ax = Axes3D(fig)
ax.plot(track[:,0], track[:,1], track[:,2],"b")

plt.show()
print("r=28,(X,Y,Z)=(0.01,0.01,0.01),t=1e4,delta_t=1e-3")

#定义x,y,z的下边界xa,ya,za，上边界xb,yb,zb
xa,ya,za=-20.0,-30.0,0.0
xb,yb,zb=20.0,30.0,50.0

#n是点的个数
n=round(T_end/T_step)

for i in range(30):
Array=[]
N=0
for i in range(n):
u,v,w=0,0,0
u=math.floor(((track[i,0])-xa)/delta)
v=math.floor(((track[i,1])-ya)/delta)
w=math.floor(((track[i,2])-za)/delta)
if [u,v,w] not in Array:
Array.append([u,v,w])
N=N+1
print(delta,N)
delta=delta/2
```

This snippet took 0.01 seconds to highlight.

Back to the Entry List or Home.