使用 Python 绘制 Poincare 截面

Question

我正要绘制下面DE的Poincare截面，这里面有一个周期势函数V(x) = - cos(x)很有意义等式.

使用时间间隔dt = 0.001的RK4计算解后，python画出如下图

但是根据教科书（J.M.T. Thompson 和 H.B. Stewart 指的是 2E），该部分看起来像

:

差别太大了。就我个人的看法，既然庞加莱截面不像作者画的那样出现，那一定是我的代码有什么错误。但是，我实际上对其他的强迫振荡DE，包括Duffing方程也做了，得到了和课本上一模一样的。所以，我想知道教科书或其他地方给出的方程式是否有错别字。我发布了我的代码，但可能很难理解。非常感谢处理它。

import numpy as np
import matplotlib.pylab as plt
import matplotlib as mpl

import sys
import time

state = [1]

def print_percent_done(index, total, state, title='Please wait'):
    percent_done2 = (index+1)/total*100
    percent_done = round(percent_done2, 1)

    print(f'\t⏳{title}: {percent_done}% done', end='\r')

    if percent_done2 > 99.9 and state[0]:
        print('\t✅'); state = [0]

####
no = 1
####

def multiple(n, q):
    m = n; i = 0
    while m >= 0:
        m -= q
        i += 1
    return min(abs(n - (i - 1)*q), abs(i*q - n))

# system(2)
#Basic info.

filename = 'sinPotentialWell'

# a = 1
# alpha = 0.01
# w = 4
w0 = .5

n = 1000000
h = .01
t_0 = 0
x_0 = 0.1
y_0 = 0
A = [(t_0, x_0, y_0)]

def f(t, x, y):
    return y
def g(t, x, y):
    return -0.5*y - np.sin(x) + 1.1*np.sin(0.5*t)

for i in range(n):
    t0 = A[i][0]; x0 = A[i][1]; y0 = A[i][2]
    
    k1 = f(t0, x0, y0)
    u1 = g(t0, x0, y0)

    k2 = f(t0 + h/2, x0 + h*k1/2, y0 + h*u1/2)
    u2 = g(t0 + h/2, x0 + h*k1/2, y0 + h*u1/2) 

    k3 = f(t0 + h/2, x0 + h*k2/2, y0 + h*u2/2)
    u3 = g(t0 + h/2, x0 + h*k2/2, y0 + h*u2/2)

    k4 = f(t0 + h, x0 + h*k3, y0 + h*u3)
    u4 = g(t0 + h, x0 + h*k3, y0 + h*u3)

    t = t0 + h
    x = x0 + (k1 + 2*k2 + 2*k3 + k4)*h/6
    y = y0 + (u1 + 2*u2 + 2*u3 + u4)*h/6
    A.append([t, x, y])
    if i%1000 == 0: print_percent_done(i, n, state, 'Solving given DE')

#phase diagram
print('showing 3d_(x, y, phi) graph')

PHI=[[]]; X=[[]]; Y=[[]]
PHI_period1 = []; X_period1 = []; Y_period1 = []

for i in range(n):
    if w0*A[i][0]%(2*np.pi) < 1 and w0*A[i-1][0]%(2*np.pi) > 6:
        PHI.append([]); X.append([]); Y.append([])
        PHI_period1.append((w0*A[i][0])%(2*np.pi)); X_period1.append(A[i][1]); Y_period1.append(A[i][2])

phi_period1 = np.array(PHI_period1); x_period1 = np.array(X_period1); y_period1 = np.array(Y_period1)

print('showing Poincare Section at phi=0')
plt.plot(x_period1, y_period1, 'gs', markersize = 2)
plt.plot()
plt.title('phi=0 Poincare Section')
plt.xlabel('x'); plt.ylabel('y')
plt.show()

Answer 1

如果把一些计算块分解出来，可以让代码更灵活，计算更直接。如果你能首先构建它，就不需要重建它。你想抓住 w0*t 是 2*pi 的倍数的点，所以只需构建时间循环，这样你就可以整合到 2*pi/w0 的块中，只记住有趣的点。

num_plot_points = 2000
h = .01
t,x,y = t_0,x_0,y_0

x_section,y_section = [],[]
T = 2*np.pi/w0
for k in range(num_plot_points):
    t = 0;
    while t < T-1.2*h:
        x,y = RK4step(t,x,y,h)
        t += h
    x,y = RK4step(t,x,y,T-t)
    if k%100 == 0: print_percent_done(k, num_plot_points, state, 'Solving given DE')
    x_section.append(x); y_section.append(y)

with RK4step 仅包含 RK4 步骤的代码。

这解不开谜团。如果您认为 x 是圆上的角度 theta（具有摩擦力的受力摆），那么面纱就会被揭开。因此，要获得具有相同空间位置的点，需要将其减少 2*pi 的倍数。这样做，

plt.plot([x%(2*np.pi) for x in x_section], y_section, 'gs', markersize = 2)

结果符合预期

使用 Python 绘制 Poincare 截面

Drawing Poincare Section using Python

python

matplotlib

differential-equations

runge-kutta

chaos