圆盘混沌的模拟

发表于 更新于

今天没事把之前参加某个比赛的部分完整了一下
顺便用刚学的求解器试了下求解微分方程，很多之前弄不出来的现在也可以弄了
第一部分是一个摆轮的混沌，第二个是倾斜摆轮后的结果

1
2
3
4
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint
from scipy.ndimage import label
1
2
3
4
5
def hundun(w,t,M,b,J,omega,g,r,k,m):
    x,v=w
    dxdt=v
    dvdt=1/J*(M*np.cos(omega*t)+m*g*r*np.sin(x)-b*v-k*x)
    return [dxdt,dvdt]
1
t=np.linspace(0,1000,10000)
1
2
3
4
5
6
7
8
M=10
b=0.1
J=5
omega=0.72
g=9.8
r=0.5
k=5
m=0.05
1
w0=(0,0)
1
result=odeint(hundun,w0,t,args=(M,b,J,omega,g,r,k,m))
1
result
array([[ 0.        ,  0.        ],
       [ 0.01006452,  0.20259081],
       [ 0.04107054,  0.42167591],
       ...,
       [-3.02715033,  0.05793392],
       [-3.02090985,  0.06682509],
       [-3.01379207,  0.07547364]])

1
result[0]
array([0., 0.])

1
2
plt.plot(result[:, 0], result[:, 1], color="blue", linewidth=0.7)
plt.show()

png

1

1
2
3
4
5
def hundun2(w,t,M,b,J,omega,g,r,k,x0,m):
    x,v=w
    dxdt=v
    dvdt=1/J*(M*np.cos(omega*t)+m*g*r*np.sin(x-x0)-b*v-k*(x))
    return [dxdt,dvdt]
1
2
3
4
5
6
7
8
9
10
11
w0=(0,0)
t=np.linspace(0,1000,100000)
M=7
b=0.5
J=5
omega=0.72
g=9.8
r=0.5
k=5
x0=0.1744
m=0.05
1
result2=odeint(hundun2,w0,t,args=(M,b,J,omega,g,r,k,x0,m))
1
2
plt.plot(result2[:, 0], result2[:, 1], color="blue", linewidth=0.7)
plt.show()

png

1

1