期末作业

未分类

在此输入正文

Lorenz系统

import matplotlib as mpl
from mpl_toolkits.mplot3d import Axes3D
 #import numpy as np
import matplotlib.pyplot as plt
 import math

xs, ys, zs = [], [], []
a=10.0
r=160
b=8.0 / 3.0

def mkPoints():
a=10.0
r=160
b=8.0 / 3.0
h=0.001
x0, y0, z0 = 0.1, 0, 0
for i in range(5000):
x1 = x0 + h a (y0 - x0)
y1 = y0 + h * (x0 * (r - z0) - y0)
z1 = z0 + h * (x0 * y0 - b* z0)
x0, y0, z0 = x1, y1, z1
xs.append(x0)
ys.append(y0)
zs.append(z0)
if name == "main":
mpl.rcParams["legend.fontsize"] = 15
fig = plt.figure()
ax = Axes3D(fig)

mkPoints()
ax.plot(xs, ys, zs, label = "Lorenz's strange attractor r=160")
ax.legend()

plt.show()

russler系统

import matplotlib as mpl
from mpl_toolkits.mplot3d import Axes3D

import numpy as np

import matplotlib.pyplot as plt
import math

xs, ys, zs = [], [], []
a=0.2
r=5.7
b=0.2

def mkPoints():
a=0.2
r=5.7
b=0.2
h=0.001
x0, y0, z0 = 0.1, 0, 0
for i in range(100000):
x1 = x0 - h * (y0 + z0)
y1 = y0 + h * (x0 +a*y0)
z1 = z0 + h * (z0*(x0-r)+b)
x0, y0, z0 = x1, y1, z1
xs.append(x0)
ys.append(y0)
zs.append(z0)
if name == "main":
mpl.rcParams["legend.fontsize"] = 15
fig = plt.figure()
ax = Axes3D(fig)

mkPoints()
ax.plot(xs, ys, zs, label = "russler's strange attractor c=5.7")
ax.legend()

plt.show()

chua系统

import matplotlib as mpl
from mpl_toolkits.mplot3d import Axes3D

import numpy as np

import matplotlib.pyplot as plt
import math

xs, ys, zs = [], [], []
p=10.0
q=14.87
m=-0.68
n=-1.27

def mkPoints():
p=10.0
q=14.87
m=-0.68
n=-1.27
h=0.001
x0, y0, z0 = 0.1, 0, 0
for i in range(100000):
x1 = x0 + h * p*(-x0+y0-m*x0-0.5*(n-m)*(abs(x0+1)-abs(x0-1)))
y1 = y0 + h * (x0-y0+z0)
z1 = z0 - h * q*y0
x0, y0, z0 = x1, y1, z1
xs.append(x0)
ys.append(y0)
zs.append(z0)
if name == "main":
mpl.rcParams["legend.fontsize"] = 15
fig = plt.figure()
ax = Axes3D(fig)

mkPoints()
ax.plot(xs, ys, zs, label = "chua's strange attractor ")
ax.legend()

plt.show()

曼德勃罗集

import numpy as np
import pylab as pl
from matplotlib import cm
def fractal(c):
z = c
for i in range(1, 1000):
if abs(z)>5: break
z = z*z+c
return i
def mandelbrot(cx, cy, d):
x0, x1, y0, y1 = cx-d, cx+d, cy-d, cy+d
y,x = np.ogrid[y0:y1:400j,x0:x1:400j]
c = x + y*1j
mandelbrot = np.frompyfunc(fractal,1,1)(c).astype(np.float)
pl.imshow(mandelbrot, cmap=cm.Greens_r, extent=[x0,x1,y0,y1])
mandelbrot(-0.5,0,1.2)
pl.show()