Как реализовать метод генерации сечений Пуанкаре для нелинейной системы ОДУ?

Question

Я пытался выяснить, как рассчитать сечения Пуанкаре для системы нелинейных ОДУ, используя документ о точной системе в качестве эталона, и боролся с 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

Чтобы получить карту Pointcaré потока ABC

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, представляют одну и ту же точку.При расчете сечения это необходимо отразить либо путем вычисления его на декартовом произведении трех кругов.Давайте останемся со вторым вариантом и выберем [-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

Чтобы получить полное представление о поперечном сечении Пуанкаре и избежать дублирования работы, используйтесетка квадратов и отметьте, если один из пересечений уже упал в нем.Новые итерации начинайте только с центров свободных квадратов.

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)

Теперь вы можете играть с плотностью сетки и длинойИнтервал интеграции, чтобы получить сюжет чуть более заполнен, но все характерные особенности уже здесь.Но я бы порекомендовал перепрограммировать это на скомпилированном языке, так как вычисления займут некоторое время.