问用欧拉法求解四维耦合系统

EN

from __future__ import division
import numpy as np    
import math    
import random   
import matplotlib.pyplot as plt    
from mpl_toolkits.mplot3d import Axes3D        
# import pdb
# pdb.set_trace()

def sys1(x, y, z, u, w , a=10, b=8.0/3.0, c=28, k1=0.4, k2=8, h=-2):

    dx = a*(y-x) + u
    dy = c*x- x*z + w    
    dz = -b*z + x*y    
    du = -h*u - x*z    
    dw = k1*x - k2*y    
    return dx, dy, dz, du, dw

t = 0.01    
itera = 2500

# Need one more for the initial values

xs = np.empty((itera+1,))    
ys = np.empty((itera+1,))    
zs = np.empty((itera+1,))    
us = np.empty((itera+1,))    
ws = np.empty((itera+1,))

# Setting initial values

xs[0], ys[0], zs[0], us[0], ws[0] = (0.1, 0.1, 0.1, 0.1, 0.1)


# Stepping through "time".

for i in range(itera):    
# Derivatives of the X, Y, Z state
    dx, dy, dz, du, dw = sys1(xs[i], ys[i], zs[i], us[i], ws[i])     

    xs[i+1] = xs[i] + (dx * t)
    ys[i+1] = ys[i] + (dy * t)
    zs[i+1] = zs[i] + (dz * t)
    us[i+1] = us[i] + (du * t)
    ws[i+1] = ws[i] + (dw * t)

fig = plt.figure()

ax = fig.gca(projection='3d')
ax.plot(xs, ys, zs)
ax.set_xlabel("X Axis")
ax.set_ylabel("Y Axis")
ax.set_zlabel("Z Axis")
ax.set_title("Lorenz Attractor")
plt.show()