CoCalc -- fmin.ipynb

Kernel: Python 2 (system-wide)

In [0]:

import numpy as np
from scipy.optimize import basinhopping
import matplotlib.pyplot as plt
import matplotlib

def rhs(X, params):
    x, y = X
    gamma, d, n, m = params
    return [gamma+d*np.sin(y-x)-np.sin(n*x), gamma+d*np.sin(x-y)-np.sin(m*y)]

Численный анализ системы $\begin{align} \dot{x} = \gamma + d \sin{(y-x)} - \sin{nx}, \\ \dot{y} = \gamma + d \sin{(x-y)} - \sin{my} \end{align}$

Построение карты режимов: цветом закодировано, есть ли в фазовом пространстве системы состояния равновесия или нет.

In [0]:

N = 100+1
grid = np.zeros([N, N])
vals = np.linspace(0.0, 5.1, N)
for i, gamma in enumerate(vals):
    doWriteOutput = True
    for j, d in enumerate(vals):
        ud = [gamma, d, 4, 5]
        if doWriteOutput:
            print ud
        rhs_eucl_sq = lambda X: np.dot(rhs(X,ud),rhs(X,ud))
        res = basinhopping(rhs_eucl_sq,[0.0,0.0],niter=500,stepsize=1.5*np.pi)
        if res.fun < 1e-12:
            if doWriteOutput:
                print "Equilibrium exists"
            grid[i][j] = 1.0
        doWriteOutput = False

%matplotlib

plt.pcolormesh(vals,vals,grid,cmap=plt.cm.get_cmap('RdBu'))
plt.xlim([min(vals), max(vals)])
plt.ylim([min(vals), max(vals)])
plt.axes().set_aspect('equal', adjustable='box')
plt.axes().set_xlabel(r'$d$, coupling strength')
plt.axes().set_ylabel(r'$\gamma$, natural frequency')
plt.axes().set_title('Regimes map: non-equilibrium motions vs equilibrium')
plt.colorbar()
plt.savefig('regimes.png')

Heatmap для нахождения координат состояния равновесия на глаз. Красный цвет соответствует минимальному значению функции (около -4), синий цвет — порядка +4. Цвет кодирует значения $\ln{(10^{-4}+P^2(x, y) + Q^2(x, y))}$ — эта функция показывает насколько мал по модулю вектор фазовой скорости системы с векторным полем $\dot{x} = P(x, y),\, \dot{y} = Q(x, y)$ .

In [0]:

# Plot equilibria

%matplotlib inline
ud = [1.2,0.4,3,4]
xs = ys = np.linspace(-np.pi,+np.pi, 2001)
res = np.zeros([len(xs),len(xs)])
for i, y in enumerate(ys):
    for j, x in enumerate(xs):
        #res[i][j] = np.dot(rhs([x,y],ud),rhs([x,y],ud))
        res[i][j] = np.log10(np.dot(rhs([x,y],ud),rhs([x,y],ud))+1e-14)
        #res[i][j] = min(np.dot(rhs([x,y],ud),rhs([x,y],ud)), 0.5)
matplotlib.rcParams['figure.figsize'] = 10, 20        
plt.pcolormesh(xs,ys,res,cmap=plt.cm.get_cmap('RdBu'))
plt.xlim([-np.pi,+np.pi])
plt.ylim([-np.pi,+np.pi])
plt.axes().set_aspect('equal', adjustable='box')
print("Min = {}".format(res.min()))
#plt.colorbar()

In [0]:

xs = np.linspace(-np.pi, +np.pi, 201)
gamma = ud[0]
matplotlib.rcParams['figure.figsize'] = 5, 10       
ys1 = gamma - np.sin(3*xs)
ys2 = gamma - np.sin(4*xs)
plt.plot(xs, ys1, color='red')
plt.plot(xs, ys2, color='green')
plt.plot([-np.pi, +np.pi], [0.,0.], 'k--')
plt.axes().set_aspect('equal', adjustable='box')

На фоне предыдущего heatmap'а изображены численно построенные фазовые траектории. Зелёные точки отмечают начальные условия, жёлтые — отдельные точки траекторий, взятые с равномерным шагом по времени, оранжевые соответствуют конечным точкам. Случайно выбирается 1000 начальных условий.

In [0]:

# Plot equilibria and sample trajectories

from scipy.integrate import odeint

def fitToStandardRange(x):
    if x < np.pi and x >= -np.pi:
        return x
    elif x >= np.pi:
        return fitToStandardRange(x-2*np.pi)
    elif x < -np.pi:
        return fitToStandardRange(x+2*np.pi)

%matplotlib inline
ud = [0.2,1,2,3]
#ud = [0.2,0.8,2,3]
#ud = [0.2,0.0,2,3]
#ud = [0.2,100.0,2,3]
xs = ys = np.linspace(-np.pi,+np.pi, 2001)
res = np.zeros([len(xs),len(xs)])
for i, y in enumerate(ys):
    for j, x in enumerate(xs):
        #res[i][j] = np.dot(rhs([x,y],ud),rhs([x,y],ud))
        res[i][j] = np.log10(np.dot(rhs([x,y],ud),rhs([x,y],ud))+1e-4)
        #res[i][j] = min(np.dot(rhs([x,y],ud),rhs([x,y],ud)), 0.5)
matplotlib.rcParams['figure.figsize'] = 20, 40        
plt.pcolormesh(xs,ys,res,cmap=plt.cm.get_cmap('RdBu'))
plt.xlim([-np.pi,+np.pi])
plt.ylim([-np.pi,+np.pi])
plt.axes().set_aspect('equal', adjustable='box')
#plt.colorbar()

n_pts = 1000
pt_xs, pt_ys = 2*np.pi*np.random.random((n_pts,))-np.pi, 2*np.pi*np.random.random((n_pts,))-np.pi
pts = zip(pt_xs, pt_ys)

rhs_vec = lambda X, t: rhs(X, ud)
time = np.linspace(0.0, 40, 1000)
for pt in pts:
    trajectory = odeint(rhs_vec,pt, time)
    trajectory = [(fitToStandardRange(x), fitToStandardRange(y)) for x, y in trajectory]
    xs, ys = zip(*trajectory)
    plt.scatter(xs,ys, color='yellow', s =0.25)
    plt.scatter([trajectory[0][0]],[trajectory[0][1]],color='green', s=25)
    plt.scatter([trajectory[-1][0]],[trajectory[-1][1]],color='orange', s=25)

In [0]: