2

我正要绘制以下 DE 的 Poincare 截面，这对于在这个方程中具有周期性势函数V(x) = - cos(x)非常有意义。

DE 我一直在尝试画庞加莱截面

使用 RK4 计算时间间隔 dt = 0.001 的解后，python 绘制的解如下图所示。

在此处输入图像描述

但根据教科书（JMT Thompson 和 HB Stewart 称为 2E），该部分看起来像

教科书中绘制的庞加莱截面：

它有很大的不同。就我个人而言，由于 Poincare 部分并没有像作者绘制的那样出现，所以我的代码中一定有一些错误。然而，我实际上对其他强迫振荡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()
4

1 回答 1

3

如果你分解出一些计算块，你可以使代码更灵活，计算更直接。如果您可以首先构建它，则无需重建某些东西。您想捕捉 的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)

仅包含RK4stepRK4 步骤的代码。

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

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

结果在预期的情节

庞加莱节

于 2021-12-10T09:50:52.123 回答