Sharedfmin.ipynbOpen in CoCalc
Author: Evgenij Grines
Views : 11
In [ ]:
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)]

Численный анализ системы

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

In [ ]:
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(104+P2(x,y)+Q2(x,y))\ln{(10^{-4}+P^2(x, y) + Q^2(x, y))} — эта функция показывает насколько мал по модулю вектор фазовой скорости системы с векторным полем x˙=P(x,y),y˙=Q(x,y)\dot{x} = P(x, y),\, \dot{y} = Q(x, y).

In [ ]:
# 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 [ ]:
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 [ ]:
# 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 [ ]: