from scipy.integrate import odeint
import numpy as np
import pylab as plt
np.random.seed(2) #為了進行一致性比較，每次執行取相同隨機數
sigma=10; rho=28; beta=8/3;
g=lambda f,t: [sigma(f[1]-f[0]), rhof[0]-f[1]-f[0]f[2],
f[0]f[1]-beta*f[2]] #定義微分方程組的右端項
s01=np.random.rand(3) #初始值
t0=np.linspace(0,50,5000)
s1=odeint(g,s01,t0) #求數值解
plt.rc('font', family='SimHei')
plt.rc('axes', unicode_minus=False)
ax=plt.subplot(121, projection='3d')
plt.plot(s1[:,0],s1[:,1],s1[:,2],'r') #畫軌線
ax.set_xlabel('$x$'); ax.set_ylabel('$y$'); ax.set_zlabel('$z$')
s02=s01+0.000001
s2 = odeint(g,s02,t0) #初值變化後，再求數值解
plt.subplot(122)
plt.plot(t0,s1[:,0]-s2[:,0],'.-') #畫x(t)的差
plt.xlabel('$t$'); plt.ylabel('$x_1(t)-x_2(t)$',rotation=90)
plt.show()
2023310143007
