如何实现为非线性 ODE 系统生成 Poincaré 截面的方法？

Question

我一直在尝试研究如何计算非线性 ODE 系统的庞加莱截面，使用关于精确系统的论文作为参考，并且一直在与 numpy 进行角力以尝试制作它运行更好。这是为了运行在有界域内。

目前我有以下代码

import numpy as np
from scipy.integrate import odeint
X = 0
Y = 1
Z = 2

def generate_poincare_map(function, initial, plane, iterations, delta):
    intersections = []
    p_i = odeint(function, initial.flatten(), [0, delta])[-1]
    for i in range(1, iterations):
        p_f = odeint(function, p_i, [i * delta, (i+1) * delta])[-1]
        if (p_f[Z] > plane) and (p_i[Z] < plane):
            intersections.append(p_i[:2])
        if (p_f[Z] > plane) and (p_i[Z] < plane):
            intersections.append(p_i[:2])
        p_i = p_f
    return np.stack(intersections)

这非常浪费，因为仅在连续的时间步长之间进行积分，而且似乎会产生不正确的结果。原始参考包括沿线的部分

而我的结果往往类似于

对于如何使这个更正确，也许更快一点，您有什么建议吗？

Answer 1

获取 ABC 流程的 Pointcaré 地图

def ABC_ode(u,t):
    A, B, C = 0.75, 1, 1 # matlab parameters
    x, y, z = u
    return np.array([
        A*np.sin(z)+C*np.cos(y), 
        B*np.sin(x)+A*np.cos(z), 
        C*np.sin(y)+B*np.cos(x)
    ])

def mysolver(u0, tspan): return odeint(ABC_ode, u0, tspan, atol=1e-10, rtol=1e-11)

你首先要明白动力系统实际上是关于单位圆上的点(cos(x),sin(x))等。因此，相差 2*pi 倍数的值表示同一点。在本节的计算中，必须通过对 3 个圆的笛卡尔积进行计算来反映这一点。让我们继续使用第二个变体，并选择 [-pi,pi] 作为基本周期，使零位置位于中心。请记住，更大的跳跃 pi 来自角度减小，而不是来自该间隔的真正交叉。

def find_crosssections(x0,y0):
    u0 = [x0,y0,0]
    px = []
    py = []
        
    u = mysolver(u0, np.arange(0, 4000, 0.5)); u0 = u[-1]
    u = np.mod(u+pi,2*pi)-pi
    x,y,z = u.T

    for k in range(len(z)-1): 
        if z[k]<=0 and z[k+1]>=0 and z[k+1]-z[k]<pi:
            # find a more exact intersection location by linear interpolation
            s = -z[k]/(z[k+1]-z[k])  # 0 = z[k] + s*(z[k+1]-z[k])
            rx, ry = (1-s)*x[k]+s*x[k+1], (1-s)*y[k]+s*y[k+1]
            px.append(rx); 
            py.append(ry);
    return px,py

要全面了解 Poincare 横截面并避免重复工作，请使用正方形网格并标记其中是否有一个交叉点。只从自由方块的中心开始新的迭代。

N=20
grid = np.zeros([N,N], dtype=int)
for i in range(N):
    for j in range(N):
        if grid[i,j]>0: continue;
        x0, y0 = (2*i+1)*pi/N-pi, (2*j+1)*pi/N-pi 
        px, py = find_crosssections(x0,y0)
        for rx,ry in zip(px,py):
            m, n = int((rx+pi)*N/(2*pi)), int((ry+pi)*N/(2*pi))
            grid[m,n]=1
        
    plt.plot(px, py, '.', ms=2)

您现在可以调整网格的密度和积分区间的长度来使绘图更完整一些，但所有的特征都已经在这里了。但我建议用编译语言重新编程，因为计算需要一些时间。

如何实现为非线性 ODE 系统生成 Poincaré 截面的方法？

How to implement a method to generate Poincaré sections for a non-linear system of ODEs?

numpy

scipy

numerical-integration

ode