Построение бассейна притяжения генератора Duffing с матплотлибом - PullRequest
Новая
       131

Построение бассейна притяжения генератора Duffing с матплотлибом

0 голосов
/ 18 октября 2018

У меня есть система обыкновенных дифференциальных уравнений, в которой оказывается два аттрактора, один в (1, 0), а другой в (-1, 0).Я хотел бы нарисовать область притяжения в декартовой координате, где есть два цвета, показывающие, какой аттрактор точка в каждой точке координаты будет заканчиваться, когда время стремится к положительной бесконечности.Однако я не знаю, как построить такой график с помощью matplotlib.Вот что я получил до сих пор: 

from scipy.integrate import ode
import numpy as np
import matplotlib.pyplot as plt
from numpy.linalg import norm

"""
The system of ODE:
x' = y
y' = x - x**3 - gamma*y
"""

# The system of equation
def f(t, r, arg):
    return [r[1], r[0] - r[0] ** 3 - arg * r[1]]

# The Jacobian matrix
def jac(t, r, arg):
    return [[0, 1], [1 - 3 * r[0] ** 2, -arg]]

# r is the vector (x,y)
# Initial condition, length of time evolution, time step, parameter gamma
def solver(r0, t0, t1, dt, gamma):

    solution = ode(f, jac).set_integrator('dopri5')

    # Set the value of gamma
    solution.set_initial_value(r0, t0).set_f_params(gamma).set_jac_params(gamma)

    return solution


# The function to find the fixed point each starting point ends at
def find_fp(r0, t0, t1, dt, gamma):
    solution = solver(r0, t0, t1, dt, gamma)
    error = 0.01
    while solution.successful():
        if norm(np.array(solution.integrate(solution.t+dt)) - np.array([1, 0])) < error:
            return 1
        elif norm(np.array(solution.integrate(solution.t+dt)) - np.array([-1, 0])) < error:
            return -1

def fp(i, j, gamma):
    t0, t1, dt = 0, 10, 0.1
    return find_fp([i, j], t0, t1, dt, gamma)

Я определил несколько функций.f - это функция, определяющая систему уравнений, jac матрица Якоби системы, которая служит параметрами для решения ОДУ с использованием метода dopri5 для scipy.integrate.ode (метод Кутта-Рунге).Функция find_fp определена для возврата аттрактора, которым будет заканчиваться точка в фазовом пространстве, возвращаемое значение 1 означает, что точка будет заканчиваться в (1, 0), а -1 в (-1).0).Функции, кажется, работают хорошо до сих пор.Тем не менее, у меня нет идей, как построить область притяжения, используя то, что я сделал с модулем matplotlib.Есть ли хорошие идеи по этому поводу?

1 Ответ

0 голосов
/ 18 октября 2018

Quick'n'dirty: выберите начальные точки, близкие к стационарным точкам, и рассчитайте решение в обратном направлении в течение некоторого времени.Постройте решения и раскрасьте их в соответствии с аттрактором. 

gamma = 1.2

def solution(x,y):
    return odeint(lambda u,t: -np.array(f(t,u,gamma)), [x,y], np.arange(0,15,0.01))

for i in range(-10,11):
    for j in range(-10,11):
        sol = solution(-1+i*1e-4, 0+j*1e-4);
        x,y = sol.T; plt.plot(x,y,'.b',ms=0.5);
        sol = solution(+1+i*1e-4, 0+j*1e-4);
        x,y = sol.T; plt.plot(x,y,'.r',ms=0.5);
plt.grid(); plt.show();

Это дает изображение

phase portrait

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

...