lsa matrix changed for new active stress

no_c_lsa.py

@@ -3,12 +3,12 @@ import matplotlib.pyplot as plt
 from matplotlib.widgets import Slider
 from matplotlib.widgets import Slider
 
 
 # real part of largest eigenvalue of STABILITY MATRIX
 # real part of largest eigenvalue of STABILITY MATRIX
-def largest_real_eigval(q, K=1.0, lamda=1.5, eta = 1, k_rho=10.0, rho_0=1, D_rho=1):
-    A = np.asmatrix([[-K*q**2/(1+eta*q**2),1j*lamda*q/(1+eta*q**2)], [1j*K*rho_0*q**3/(1+eta*q**2),-D_rho*q**2 -k_rho + rho_0*lamda*q**2/(1+eta*q**2)]]) 
+def largest_real_eigval(q, K=1.0, lamda=1.5, eta = 1, k_rho=10.0, rho_0=1, D_rho=1, rho_s=1.0):
+    A = np.asmatrix([[-K*q**2/(1+eta*q**2),1j*lamda*q*rho_s/((1+eta*q**2)*(rho_0 + rho_s)**2)], [1j*K*rho_0*q**3/(1+eta*q**2),-D_rho*q**2 -k_rho + rho_0*lamda*q**2*rho_s/((1+eta*q**2)*(rho_0 + rho_s)**2)]]) 
     lamda = np.real(np.linalg.eigvals(A))
     lamda = np.real(np.linalg.eigvals(A))
     return lamda.max()
     return lamda.max()
 
 
-kyu = np.linspace(0,6*np.pi,100)
+kyu = np.linspace(0,2*np.pi,1000)
 
 
 fig, ax = plt.subplots(1, figsize=(8,6))
 fig, ax = plt.subplots(1, figsize=(8,6))
 fig.subplots_adjust(left=0.15, bottom=0.3, right=0.98, top=0.95,
 fig.subplots_adjust(left=0.15, bottom=0.3, right=0.98, top=0.95,
@@ -18,13 +18,14 @@ ax.set_ylabel(r'$Re[\, \lambda(q) \, ]_{\rm max}$')
 ax.axhline(y=0, color='black')
 ax.axhline(y=0, color='black')
 
 
 D_rho_0 = 1.0
 D_rho_0 = 1.0
-lamda0 = 6.0
+lamda0 = 20.0
 K0 = 1.0
 K0 = 1.0
 k_rho0 = 1.0
 k_rho0 = 1.0
 eta0 = 1.0
 eta0 = 1.0
 rho_00 = 1.0
 rho_00 = 1.0
+rho_s0 = 1.0
 
 
-lamda = np.array([largest_real_eigval(q, K=K0, lamda=lamda0, eta=eta0, k_rho=k_rho0, rho_0=rho_00, D_rho=D_rho_0) for q in kyu])
+lamda = np.array([largest_real_eigval(q, K=K0, lamda=lamda0, eta=eta0, k_rho=k_rho0, rho_0=rho_00, D_rho=D_rho_0, rho_s = rho_s0) for q in kyu])
 lambda_plot, = ax.plot(kyu, lamda)
 lambda_plot, = ax.plot(kyu, lamda)
 lamda_min, lamda_max = min(0, lamda.min()), max(0, lamda.max())
 lamda_min, lamda_max = min(0, lamda.min()), max(0, lamda.max())
 ax.set_ylim(lamda_min, lamda_max)
 ax.set_ylim(lamda_min, lamda_max)
@@ -36,18 +37,20 @@ ax_eta = plt.axes([0.4, 0.15, 0.2, 0.02])
 ax_k_rho = plt.axes([0.4, 0.20, 0.2, 0.02])
 ax_k_rho = plt.axes([0.4, 0.20, 0.2, 0.02])
 ax_rho_0  = plt.axes([0.7, 0.15, 0.2, 0.02])
 ax_rho_0  = plt.axes([0.7, 0.15, 0.2, 0.02])
 ax_D_rho  = plt.axes([0.7, 0.20, 0.2, 0.02])
 ax_D_rho  = plt.axes([0.7, 0.20, 0.2, 0.02])
+ax_rho_s  = plt.axes([0.7, 0.10, 0.2, 0.02])
 
 
 # sliders for controlling parameters
 # sliders for controlling parameters
 s_K  = Slider(ax_K, r'$K$', valmin=0.0, valmax=10.0, valinit=K0, valstep=0.001)
 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_lamda  = Slider(ax_lamda, r'$\lambda$', valmin=0.0, valmax=100.0, valinit=lamda0, valstep=0.001)
 s_k_rho  = Slider(ax_k_rho, r'$k_{\rho}$', valmin=0.0, valmax=10.0, valinit=k_rho0, valstep=0.001)
 s_k_rho  = Slider(ax_k_rho, r'$k_{\rho}$', valmin=0.0, valmax=10.0, valinit=k_rho0, valstep=0.001)
 s_rho_0  = Slider(ax_rho_0, r'$\rho_0$', valmin=0.0, valmax=10.0, valinit=rho_00, valstep=0.001)
 s_rho_0  = Slider(ax_rho_0, r'$\rho_0$', valmin=0.0, valmax=10.0, valinit=rho_00, valstep=0.001)
 s_eta  = Slider(ax_eta, r'$\eta$', valmin=0.0, valmax=10.0, valinit=eta0, valstep=0.001)
 s_eta  = Slider(ax_eta, r'$\eta$', valmin=0.0, valmax=10.0, valinit=eta0, valstep=0.001)
 s_D_rho  = Slider(ax_D_rho, r'$D_{\rho}$', valmin=0.0, valmax=10.0, valinit=D_rho_0, valstep=0.001)
 s_D_rho  = Slider(ax_D_rho, r'$D_{\rho}$', valmin=0.0, valmax=10.0, valinit=D_rho_0, valstep=0.001)
+s_rho_s  = Slider(ax_rho_s, r'$\rho_s$', valmin=0.0, valmax=10.0, valinit=rho_s0, valstep=0.001)
 
 
 # slider update function
 # slider update function
 def update(val):
 def update(val):
-    lamda = np.array([largest_real_eigval(q, K=s_K.val, lamda=s_lamda.val, eta=s_eta.val, k_rho=s_eta.val, rho_0=s_rho_0.val, D_rho=s_D_rho.val)
+    lamda = np.array([largest_real_eigval(q, K=s_K.val, lamda=s_lamda.val, eta=s_eta.val, k_rho=s_eta.val, rho_0=s_rho_0.val, D_rho=s_D_rho.val, rho_s = s_rho_s.val)
                             for q in kyu])
                             for q in kyu])
     lambda_plot.set_ydata(lamda)
     lambda_plot.set_ydata(lamda)
     lamda_min, lamda_max = min(0, lamda.min()), max(0, lamda.max())
     lamda_min, lamda_max = min(0, lamda.min()), max(0, lamda.max())
@@ -60,5 +63,6 @@ s_eta.on_changed(update)
 s_k_rho.on_changed(update)
 s_k_rho.on_changed(update)
 s_rho_0.on_changed(update)
 s_rho_0.on_changed(update)
 s_D_rho.on_changed(update)
 s_D_rho.on_changed(update)
+s_rho_s.on_changed(update)
 
 
 plt.show()
 plt.show()