with signaling -- linear stability

with_c_lsa.py

+import numpy as np
+import matplotlib.pyplot as plt
+from matplotlib.widgets import Slider
+
+# real part of largest eigenvalue of STABILITY MATRIX
+def largest_real_eigval(q, K=1.0, lamda=1.5, eta = 1, krho=10.0, rho0 = 1, c0 = 1.0, kc = 1.0, Drho = 1.0, Dc = 1.0):
+    A = np.asmatrix([[-K * q**2/(1 + eta * q**2), 1j * c0 * lamda/(1 + eta * q**2), 1j * rho0 * lamda/(1 + eta * q**2)],
+                    [1j * q**3 * rho0 * K/(1 + eta * q**2), -Drho * q**2 -krho + q**2 * c0 * rho0 * lamda/(1 + eta * q**2), q**2 * rho0**2 * lamda/(1 + eta * q**2)],
+                    [1j * q**3 * c0 * K/(1 + eta * q**2), q**2 * c0**2 * lamda/(1 + eta * q**2), -kc + q**2 * c0 * rho0 * lamda/(1 + eta * q**2) - q**2 * Dc]]) 
+    lamda = np.real(np.linalg.eigvals(A))
+    return lamda.max()
+
+kyu = np.linspace(0,100,1000)
+
+fig, ax = plt.subplots(1, figsize=(8,6))
+fig.subplots_adjust(left=0.15, bottom=0.3, right=0.98, top=0.95,
+        wspace=0.1, hspace=0.2)
+ax.set_xlabel(r'$q$')
+ax.set_ylabel(r'$Re[\, \lambda(q) \, ]_{\rm max}$')
+ax.axhline(y=0, color='black')
+
+lamda0 = 3.88
+K0 = 0.94
+krho0 = 2.81
+eta0 = 3.12
+rho00 = 1.0
+Dc0 = 1.12
+Drho0 = 0.0
+kc0 = 1.38
+c00 = 2.12
+
+lamda = np.array([largest_real_eigval(q, K = K0, lamda = lamda0, eta = eta0, krho = krho0, rho0 = rho00, c0 = c00, kc = kc0, Drho = Drho0, Dc = Dc0) for q in kyu])
+lambda_plot, = ax.plot(kyu, lamda)
+lamda_min, lamda_max = min(0, lamda.min()), max(0, lamda.max())
+ax.set_ylim(lamda_min, lamda_max)
+
+# axes for controlling parameter sliders
+ax_K  = plt.axes([0.1, 0.15, 0.2, 0.02])
+ax_lamda  = plt.axes([0.1, 0.20, 0.2, 0.02])
+ax_eta = plt.axes([0.4, 0.15, 0.2, 0.02])
+ax_krho = plt.axes([0.4, 0.20, 0.2, 0.02])
+ax_rho0  = plt.axes([0.7, 0.15, 0.2, 0.02])
+ax_Dc  = plt.axes([0.7, 0.20, 0.2, 0.02])
+ax_kc = plt.axes([0.1, 0.10, 0.2, 0.02])
+ax_c0 = plt.axes([0.4, 0.10, 0.2, 0.02])
+ax_Drho = plt.axes([0.7, 0.10, 0.2, 0.02])
+
+# sliders for controlling parameters
+s_K  = Slider(ax_K, r'$K$', valmin = 0.0, valmax = 10.0, valinit = K0, valstep = 0.001)
+s_lamda  = Slider(ax_lamda, r'$\lambda$', valmin = 0.0, valmax = 20.0, valinit = lamda0, valstep = 0.001)
+s_eta  = Slider(ax_eta, r'$\eta$', valmin = 0.0, valmax = 10.0, valinit = eta0, valstep = 0.001)
+s_krho  = Slider(ax_krho, r'$k_{\rho}$', valmin = 0.0, valmax = 10.0, valinit = krho0, valstep = 0.001)
+s_rho0  = Slider(ax_rho0, r'$\rho_0$', valmin = 0.0, valmax = 10.0, valinit = rho00, valstep = 0.001)
+s_Dc  = Slider(ax_Dc, r'$D_c$', valmin = 0.0, valmax = 10.0, valinit = Dc0, valstep = 0.001)
+s_kc  = Slider(ax_kc, r'$k_c$', valmin = 0.0, valmax = 10.0, valinit = kc0, valstep = 0.001)
+s_c0  = Slider(ax_c0, r'$c_0$', valmin = 0.0, valmax = 10.0, valinit = c00, valstep = 0.001)
+s_Drho  = Slider(ax_Drho, r'$D_\rho{}$', valmin = 0.0, valmax = 10.0, valinit = Drho0, valstep = 0.001)
+
+# slider update function
+def update(val):
+    lamda = np.array([largest_real_eigval(q, K=s_K.val, lamda = s_lamda.val, eta = s_eta.val, krho = s_krho.val, rho0 = s_rho0.val, c0 = s_c0.val, kc = s_kc.val, Drho = s_Drho.val, Dc = s_Dc.val)
+                            for q in kyu])
+    lambda_plot.set_ydata(lamda)
+    lamda_min, lamda_max = min(0, lamda.min()), max(0, lamda.max())
+    ax.set_ylim(lamda_min, lamda_max)
+    fig.canvas.draw()
+
+s_K.on_changed(update)
+s_lamda.on_changed(update)
+s_eta.on_changed(update)
+s_krho.on_changed(update)
+s_rho0.on_changed(update)
+s_Dc.on_changed(update)
+s_Drho.on_changed(update)
+s_kc.on_changed(update)
+s_c0.on_changed(update)
+
+plt.show()