Drawing Poincare Section using Python-CodePudding

I was about to plot a Poincare section of the following DE, which is quite meaningful to have a periodic potential function V(x) = - cos(x) in this equation.

After calculating the solution using RK4 with time interval dt = 0.001, the one that python drew was as the following plot.

But according to the textbook(referred to 2E by J.M.T. Thompson and H.B. Stewart), the section would look like as

it has so much difference. For my personal opinion, since Poincare section does not appear as what writers draw, there must be some error in my code. However, I actually done for other forced oscillation DE, including Duffing's equation, and obtained the identical one as those in the textbook. So, I was wodering if there are some typos in the equation given by the textbook, or somewhere else. I posted my code, but might be quite messy to understand. So appreicate dealing with it.

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 i00 == 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()

CodePudding user response：

If you factor out some of the computation blocks, you can make the code more flexible and computations more direct. No need to reconstruct something if you can construct it in the first place. You want to catch the points where w0*t is a multiple of 2*pi, so just construct the time loops so you integrate in chunks of 2*pi/w0 and only remember the interesting points.

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 k0 == 0: print_percent_done(k, num_plot_points, state, 'Solving given DE')
    x_section.append(x); y_section.append(y)

with RK4step just containing the code of the RK4 step.

This will not solve the mystery. The veil gets lifted if you consider that x is the angle theta (of a forced pendulum with friction) on a circle. Thus to get points with the same spacial location it needs to be reduced by multiples of 2*pi. Doing that,

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

results in the expected plot