CoCalc Shared Filesgrim / src / problem / linear_modes / Balbusaur.sagews
Authors: Mani Chandra, Francois Foucart
Views : 7
%auto
typeset_mode(True)


## $\mathtt{balbusaur}$

#### Initialize the operators $\mathtt{d\_dt()}$, $\mathtt{d\_dX1()}$ and $\mathtt{d\_dX2()}$ acting on variables:

1. Density $\rho$
2. Internal energy $u$
3. Velocity in $X^1$ direction $u^1$
4. Velocity in $X^2$ direction $u^2$
5. Velocity in $X^3$ direction $u^3$
6. Magnetic field in $X^1$ direction $B^1$
7. Magnetic field in $X^2$ direction $B^2$
8. Magnetic field in $X^3$ direction $B^3$
9. Heat flux magnitude $q$
10. Pressure anisotropy magnitude $dp$
• Mean variables : $\rho_0$, $u_0$, $u^1_0$, $u^2_0$, $u^3_0$, $B^1_0$, $B^2_0$, $B^3_0$, $q_0$, $dp_0$
• Perturbed variables : $\delta_\rho$, $\delta_u$, $\delta_{u^1}$, $\delta_{u^2}$, $\delta_{u^3}$, $\delta_{B^1}$, $\delta_{B^2}$, $\delta_{B^3}$, $\delta_q$, $\delta_{dp}$
def linearize(term):
return taylor(term, (delta_rho, 0), \
(delta_u, 0),   \
(delta_u1, 0),  \
(delta_u2, 0),  \
(delta_u3, 0),  \
(delta_B1, 0),  \
(delta_B2, 0),  \
(delta_B3, 0),  \
(delta_q, 0), \
(delta_dp, 0), \
(delta_rho_dt, 0), \
(delta_u_dt, 0),   \
(delta_u1_dt, 0),  \
(delta_u2_dt, 0),  \
(delta_u3_dt, 0),  \
(delta_B1_dt, 0),  \
(delta_B2_dt, 0),  \
(delta_B3_dt, 0),  \
(delta_q_dt, 0), \
(delta_dp_dt, 0), 1 \
).simplify_full()

def d_dX1(term):
term  = Expression(SR, linearize(term))

expr =   term.coefficient(delta_rho) * I * k1 * delta_rho \
+ term.coefficient(delta_u)   * I * k1 * delta_u   \
+ term.coefficient(delta_u1)  * I * k1 * delta_u1  \
+ term.coefficient(delta_u2)  * I * k1 * delta_u2  \
+ term.coefficient(delta_u3)  * I * k1 * delta_u3  \
+ term.coefficient(delta_B1)  * I * k1 * delta_B1  \
+ term.coefficient(delta_B2)  * I * k1 * delta_B2  \
+ term.coefficient(delta_B3)  * I * k1 * delta_B3  \
+ term.coefficient(delta_q)   * I * k1 * delta_q \
+ term.coefficient(delta_dp)  * I * k1 * delta_dp

return expr

def d_dX2(term):
term  = Expression(SR, linearize(term))

expr =   term.coefficient(delta_rho) * I * k2 * delta_rho \
+ term.coefficient(delta_u)   * I * k2 * delta_u   \
+ term.coefficient(delta_u1)  * I * k2 * delta_u1  \
+ term.coefficient(delta_u2)  * I * k2 * delta_u2  \
+ term.coefficient(delta_u3)  * I * k2 * delta_u3  \
+ term.coefficient(delta_B1)  * I * k2 * delta_B1  \
+ term.coefficient(delta_B2)  * I * k2 * delta_B2  \
+ term.coefficient(delta_B3)  * I * k2 * delta_B3  \
+ term.coefficient(delta_q)   * I * k2 * delta_q \
+ term.coefficient(delta_dp)  * I * k2 * delta_dp

return expr

def d_dt(term):
term  = Expression(SR, linearize(term))

expr =   term.coefficient(delta_rho) * delta_rho_dt \
+ term.coefficient(delta_u)   * delta_u_dt   \
+ term.coefficient(delta_u1)  * delta_u1_dt  \
+ term.coefficient(delta_u2)  * delta_u2_dt  \
+ term.coefficient(delta_u3)  * delta_u3_dt  \
+ term.coefficient(delta_B1)  * delta_B1_dt  \
+ term.coefficient(delta_B2)  * delta_B2_dt  \
+ term.coefficient(delta_B3)  * delta_B3_dt  \
+ term.coefficient(delta_q)   * delta_q_dt \
+ term.coefficient(delta_dp)  * delta_dp_dt

return expr



#### Options:

1. $\mathtt{EVOLVE\_B\_FIELDS}$ : 0 or 1
2. $\mathtt{CONDUCTION}$ : 0 or 1
3. $\mathtt{VISCOSITY}$ : 0 or 1
4. $\mathtt{FAKE\_EMHD}$ : 0 or 1
5. $\mathtt{TURN\_OFF\_MEAN\_B2}$ : 0 or 1
6. $\mathtt{TURN\_OFF\_K2\_PERTURBATIONS}$ : 0 or 1
7. $\mathtt{PRINCIPAL\_COEFFICIENTS}$ : 0 or 1
# Spatiotemporal variables
t, omega, k1, k2 = var('t, omega, k1, k2')

# Constants:
# kappa : Heat conductivity
# eta   : shear viscosity
# tau   : relaxation time scale
# phi   : dimensionless coefficient for conduction
# psi   : dimensionless coefficient for viscosity
Gamma, kappa, eta, tau, phi, psi = var('Gamma, kappa, eta, tau, phi, psi')

# Background mean values: Symbolic variables
rho0, u0, u10, u20, u30, B10, B20, B30, q0, dp0 = var('rho0, u0, u10, u20, u30, B10, B20, B30, q0, dp0')

# Perturbations in space
delta_rho, delta_u, delta_u1, delta_u2, delta_u3, delta_B1, delta_B2, delta_B3, delta_q, delta_dp = \
var('delta_rho, delta_u, delta_u1, delta_u2, delta_u3, delta_B1, delta_B2, delta_B3, delta_q, delta_dp')

# Perturbations in time
delta_rho_dt, delta_u_dt, delta_u1_dt, delta_u2_dt, delta_u3_dt, delta_B1_dt, delta_B2_dt, delta_B3_dt, delta_q_dt, delta_dp_dt = \
var('delta_rho_dt, delta_u_dt, delta_u1_dt, delta_u2_dt, delta_u3_dt, delta_B1_dt, delta_B2_dt, delta_B3_dt, delta_q_dt, delta_dp_dt')

# Inputs:

FAKE_EMHD                 = 0
EVOLVE_B_FIELDS           = 0
CONDUCTION                = 1
VISCOSITY                 = 0
PRINCIPAL_COEFFICIENTS    = 0
TURN_ON_MEAN_B2           = 0
TURN_ON_K2_PERTURBATIONS  = 0
TURN_ON_U2                = 0
TURN_ON_U3                = 0

# Equilibrium states:
# rho0 and u0 are NEVER set to zero
u10 = 0
u20 = 0
u30 = 0
B30 = 0
#q0  = 0
dp0 = 0
if (TURN_ON_MEAN_B2==0):
B20 = 0
if (TURN_ON_K2_PERTURBATIONS==0):
k2 = 0

rho = rho0 + delta_rho
u   = u0   + delta_u
u1  = u10  + delta_u1
u2  = u20  + delta_u2
u3  = u30  + delta_u3
B1  = B10  + delta_B1
B2  = B20  + delta_B2
B3  = B30  + delta_B3
q   = q0   + CONDUCTION*delta_q
dp  = dp0  + VISCOSITY*delta_dp

# Introducing:
# chi : Thermal diffusivity
# nu  : kinematic viscosity
# cs  : sound speed
P   = (Gamma - 1)*u
T   = P/rho
cs  = sqrt(Gamma * P/ (rho + Gamma*u) )
chi   = phi * cs**2 * tau
nu    = psi * cs**2 * tau
kappa = rho * chi
eta   = rho * nu

# Inputs for numerical diagonalization for finite k modes
k1_num    = 1
k2_num    = 0*pi
rho0_num = 1
u0_num   = 0.1
u10_num  = 0
u20_num  = 0
u30_num  = 0
B10_num  = 0.1
B20_num  = 0.
B30_num  = 0
q0_num   = 0.01
dp0_num  = 0

Gamma_num = 5/3
tau_num   = 1
phi_num   = 1
psi_num   = 1


#### All the physics is below

gcon = Matrix([ [-1, 0, 0, 0],
[0, 1, 0, 0],
[0, 0, 1, 0],
[0, 0, 0, 1]
]
)

gcov = gcon.inverse()

gamma = sqrt(1 +    gcov[1][1]*u1*u1 + gcov[2][2]*u2*u2 + gcov[3][3]*u3*u3
+ 2*(gcov[1][2]*u1*u2 + gcov[1][3]*u1*u3 + gcov[2][3]*u2*u3)
)

ucon = [gamma, u1, u2, u3]
ucov = [-gamma, u1, u2, u3]

bcon0 = B1*ucov[1] + B2*ucov[2] + B3*ucov[3]
bcon1 = (B1 + bcon0*ucon[1])/ucon[0]
bcon2 = (B2 + bcon0*ucon[2])/ucon[0]
bcon3 = (B3 + bcon0*ucon[3])/ucon[0]

bcon = [bcon0, bcon1, bcon2, bcon3]
bcov = [-bcon0, bcon1, bcon2, bcon3]

bsqr = bcon[0]*bcov[0] + bcon[1]*bcov[1] + bcon[2]*bcov[2] + bcon[3]*bcov[3]

def delta(mu, nu):
if (mu==nu):
return 1
else:
return 0

def TUpDown(mu, nu):

return   (rho + u + P + bsqr + FAKE_EMHD*bsqr/6)*ucon[mu]*ucov[nu] \
+ (P + bsqr/2 + FAKE_EMHD*bsqr/6)*delta(mu, nu) - (1 + FAKE_EMHD/2)*bcon[mu]*bcov[nu] \
+ CONDUCTION*q/sqrt(bsqr) * (bcon[mu]*ucov[nu] + ucon[mu]*bcov[nu]) \
- VISCOSITY *dp/bsqr      * (bcon[mu]*bcov[nu]) + dp/3*(ucon[mu]*ucov[nu] + delta(mu, nu))

def acon(mu):
return linearize(ucon[0]*d_dt(ucon[mu]) + ucon[1]*d_dX1(ucon[mu]) + ucon[2]*d_dX2(ucon[mu]))

def qconEckart(mu):
acov = [-acon(0), acon(1), acon(2), acon(3)]

ans = -kappa*(ucon[mu]*ucon[0] + gcon[mu, 0])*(d_dt(T) +  T*acov[0]) \
-kappa*(ucon[mu]*ucon[1] + gcon[mu, 1])*(d_dX1(T) + T*acov[1]) \
-kappa*(ucon[mu]*ucon[2] + gcon[mu, 2])*(d_dX2(T) + T*acov[2])

return linearize(ans)

Eqn_rho = linearize(d_dt(rho*ucon[0])   + d_dX1(rho*ucon[1])   + d_dX2(rho*ucon[2]))
Eqn_u   = linearize(d_dt(TUpDown(0, 0)) + d_dX1(TUpDown(1, 0)) + d_dX2(TUpDown(2, 0)))
Eqn_u1  = linearize(d_dt(TUpDown(0, 1)) + d_dX1(TUpDown(1, 1)) + d_dX2(TUpDown(2, 1)))
Eqn_u2  = linearize(d_dt(TUpDown(0, 2)) + d_dX1(TUpDown(1, 2)) + d_dX2(TUpDown(2, 2)))
Eqn_u3  = linearize(d_dt(TUpDown(0, 3)) + d_dX1(TUpDown(1, 3)) + d_dX2(TUpDown(2, 3)))

Eqn_B1  = linearize(d_dt(B1) + d_dX2(bcon[1]*ucon[2] - bcon[2]*ucon[1]) )
Eqn_B2  = linearize(d_dt(B2) + d_dX1(bcon[2]*ucon[1] - bcon[1]*ucon[2]) )
Eqn_B3  = linearize(d_dt(B3) + d_dX1(bcon[3]*ucon[1] - bcon[1]*ucon[3]) + d_dX2(bcon[3]*ucon[2] - bcon[2]*ucon[3]) )

q_relaxed  = (bcov[0]*qconEckart(0) + bcov[1]*qconEckart(1) + bcov[2]*qconEckart(2) + bcov[3]*qconEckart(3) )/sqrt(bsqr)
dp_relaxed = 0
for nu in xrange(4):
dp_relaxed = dp_relaxed + 3*eta/bsqr * (bcon[nu]* (bcon[0]*d_dt(ucov[nu]) + bcon[1]*d_dX1(ucov[nu]) + bcon[2]*d_dX2(ucov[nu])) )

dp_relaxed = dp_relaxed - eta*(d_dt(ucon[0]) + d_dX1(ucon[1]) + d_dX2(ucon[2]) )

beta1  = tau/(kappa*T)
Eqn_q  = linearize(ucon[0]*d_dt(q)  + ucon[1]*d_dX1(q)  + ucon[2]*d_dX2(q)  + (PRINCIPAL_COEFFICIENTS*q  - q_relaxed)/tau
+ (q*T/(2*beta1))*(d_dt(beta1*ucon[0]/T) + d_dX1(beta1*ucon[1]/T) + d_dX2(beta1*ucon[2]/T))
)
Eqn_dp = linearize(ucon[0]*d_dt(dp) + ucon[1]*d_dX1(dp) + ucon[2]*d_dX2(dp) + (PRINCIPAL_COEFFICIENTS*dp - dp_relaxed)/tau)

Eqns          = [Eqn_rho==0, Eqn_u==0, Eqn_u1==0]
delta_vars    = [delta_rho, delta_u, delta_u1]
delta_vars_dt = [delta_rho_dt, delta_u_dt, delta_u1_dt]

if (TURN_ON_U2):
Eqns          = Eqns          + [Eqn_u2==0]
delta_vars    = delta_vars    + [delta_u2]
delta_vars_dt = delta_vars_dt + [delta_u2_dt]

if (TURN_ON_U3):
Eqns          = Eqns          + [Eqn_u3==0]
delta_vars    = delta_vars    + [delta_u3]
delta_vars_dt = delta_vars_dt + [delta_u3_dt]

if (EVOLVE_B_FIELDS):
Eqns          = Eqns          + [Eqn_B1==0,   Eqn_B2==0,   Eqn_B3==0]
delta_vars    = delta_vars    + [delta_B1,    delta_B2,    delta_B3]
delta_vars_dt = delta_vars_dt + [delta_B1_dt, delta_B2_dt, delta_B3_dt]

if (CONDUCTION):
Eqns          = Eqns          + [Eqn_q==0]
delta_vars    = delta_vars    + [delta_q]
delta_vars_dt = delta_vars_dt + [delta_q_dt]

if (VISCOSITY):
Eqns          = Eqns          + [Eqn_dp==0]
delta_vars    = delta_vars    + [delta_dp]
delta_vars_dt = delta_vars_dt + [delta_dp_dt]

solutions = solve(Eqns, delta_vars_dt, solution_dict=True)

solns_delta_vars_dt = []
for dvar_dt in delta_vars_dt:
solns_delta_vars_dt.append(solutions[0][dvar_dt])

M = jacobian(solns_delta_vars_dt, delta_vars)
M = M.apply_map(lambda x : x.simplify_full())

pretty_print("Linearized system : ", )
print("\n")
pretty_print(Matrix(delta_vars_dt).transpose(), " = ", M, Matrix(delta_vars).transpose())
print("\n\n")
pretty_print("Analytic eigenvalues and eigenvectors in the $k_1, k_2 \\rightarrow 0$ limit : ", )
M.subs(k1=0, k2=0).eigenvectors_right()

# Numerical diagonalization:

M_numerical = M.subs(rho0=rho0_num, u0=u0_num, u10=u10_num, u20=u20_num, u30=u30_num, \
B10=B10_num, B20=B20_num, B30=B30_num, q0=q0_num,   dp0=dp0_num, \
Gamma=Gamma_num, tau=tau_num, phi=phi_num, psi=psi_num, \
k1=k1_num, k2=k2_num \
)

M_numerical = M_numerical.change_ring(CDF)
eigenvecs   = M_numerical.eigenvectors_right()
#
rho_index = 0
u_index   = 1
u1_index  = 2
#u2_index  = 3
#u3_index  = 4
if (CONDUCTION==1 and VISCOSITY==0):
q_index = 3
if (CONDUCTION==0 and VISCOSITY==1):
dp_index = 5
if (CONDUCTION==1 and VISCOSITY==1):
q_index = 5
dp_index = 6
#
if (EVOLVE_B_FIELDS):
b1_index  = 5
b2_index  = 6
b3_index  = 7
#
if (CONDUCTION==1 and VISCOSITY==0):
q_index = 8

if (CONDUCTION==0 and VISCOSITY==1):
dp_index = 8

if (CONDUCTION==1 and VISCOSITY==1):
q_index  = 8
dp_index = 9

pretty_print("Numerical eigenvalues and eigenvectors for $k_1 =$", k1_num, " , $k_2$ = ", k2_num, ":\n")

for i in xrange(len(eigenvecs)):
print("--------------------------")
print("Eigenvalue   = ", eigenvecs[i][0])
print(delta_rho,  " = ", eigenvecs[i][1][0][0])
print(delta_u,    " = ", eigenvecs[i][1][0][1])
print(delta_u1,   " = ", eigenvecs[i][1][0][2])
#print(delta_u2,   " = ", eigenvecs[i][1][0][3])
#print(delta_u3,   " = ", eigenvecs[i][1][0][4])

if (EVOLVE_B_FIELDS):
print(delta_B1, " = ", eigenvecs[i][1][0][5])
print(delta_B2, " = ", eigenvecs[i][1][0][6])
print(delta_B3, " = ", eigenvecs[i][1][0][7])

if (CONDUCTION):
print(delta_q, " = ", eigenvecs[i][1][0][q_index])

if (VISCOSITY):
print(delta_dp," = ", eigenvecs[i][1][0][dp_index])

Linearized system :
$\displaystyle \left(\begin{array}{r} \delta_{\rho_{\mathit{dt}}} \\ \delta_{u_{\mathit{dt}}} \\ \delta_{\mathit{u1}_{\mathit{dt}}} \\ \delta_{q_{\mathit{dt}}} \end{array}\right)$ = $\displaystyle \left(\begin{array}{rrrr} 0 & 0 & -i \, k_{1} \rho_{0} & 0 \\ \frac{{\left(2 i \, \Gamma^{3} - 4 i \, \Gamma^{2} + 2 i \, \Gamma\right)} k_{1} \phi q_{0} u_{0}^{3}}{2 \, \Gamma \rho_{0}^{2} u_{0}^{2} + {\left(\Gamma^{2} - {\left(\Gamma^{3} - 2 \, \Gamma^{2} + \Gamma\right)} \phi\right)} \rho_{0} u_{0}^{3} - 3 \, q_{0}^{2} \rho_{0}^{2} - {\left(2 \, \Gamma q_{0}^{2} \rho_{0} - \rho_{0}^{3}\right)} u_{0}} & -\frac{{\left(-2 i \, \Gamma + 2 i\right)} k_{1} q_{0} \rho_{0} u_{0} + {\left({\left(2 i \, \Gamma^{3} - 4 i \, \Gamma^{2} + 2 i \, \Gamma\right)} k_{1} \phi + {\left(-2 i \, \Gamma^{2} + 2 i \, \Gamma\right)} k_{1}\right)} q_{0} u_{0}^{2}}{2 \, \Gamma \rho_{0} u_{0}^{2} + {\left(\Gamma^{2} - {\left(\Gamma^{3} - 2 \, \Gamma^{2} + \Gamma\right)} \phi\right)} u_{0}^{3} - 3 \, q_{0}^{2} \rho_{0} - {\left(2 \, \Gamma q_{0}^{2} - \rho_{0}^{2}\right)} u_{0}} & \frac{-2 i \, \Gamma^{2} k_{1} \rho_{0} u_{0}^{3} + 5 i \, k_{1} q_{0}^{2} \rho_{0} u_{0} - {\left(i \, \Gamma^{3} k_{1} - {\left(i \, \Gamma^{4} - 2 i \, \Gamma^{3} + i \, \Gamma^{2}\right)} k_{1} \phi\right)} u_{0}^{4} + {\left(4 i \, \Gamma k_{1} q_{0}^{2} - i \, \Gamma k_{1} \rho_{0}^{2}\right)} u_{0}^{2}}{2 \, \Gamma \rho_{0} u_{0}^{2} + {\left(\Gamma^{2} - {\left(\Gamma^{3} - 2 \, \Gamma^{2} + \Gamma\right)} \phi\right)} u_{0}^{3} - 3 \, q_{0}^{2} \rho_{0} - {\left(2 \, \Gamma q_{0}^{2} - \rho_{0}^{2}\right)} u_{0}} & \frac{-2 i \, \Gamma k_{1} \rho_{0} u_{0}^{2} - i \, k_{1} \rho_{0}^{2} u_{0} - {\left(i \, \Gamma^{2} k_{1} - {\left(i \, \Gamma^{3} - 2 i \, \Gamma^{2} + i \, \Gamma\right)} k_{1} \phi\right)} u_{0}^{3}}{2 \, \Gamma \rho_{0} u_{0}^{2} + {\left(\Gamma^{2} - {\left(\Gamma^{3} - 2 \, \Gamma^{2} + \Gamma\right)} \phi\right)} u_{0}^{3} - 3 \, q_{0}^{2} \rho_{0} - {\left(2 \, \Gamma q_{0}^{2} - \rho_{0}^{2}\right)} u_{0}} \\ \frac{{\left(-i \, \Gamma^{3} + 2 i \, \Gamma^{2} - i \, \Gamma\right)} k_{1} \phi u_{0}^{3}}{2 \, \Gamma \rho_{0}^{2} u_{0}^{2} + {\left(\Gamma^{2} - {\left(\Gamma^{3} - 2 \, \Gamma^{2} + \Gamma\right)} \phi\right)} \rho_{0} u_{0}^{3} - 3 \, q_{0}^{2} \rho_{0}^{2} - {\left(2 \, \Gamma q_{0}^{2} \rho_{0} - \rho_{0}^{3}\right)} u_{0}} & \frac{{\left(-i \, \Gamma + i\right)} k_{1} \rho_{0} u_{0} + {\left({\left(i \, \Gamma^{3} - 2 i \, \Gamma^{2} + i \, \Gamma\right)} k_{1} \phi + {\left(-i \, \Gamma^{2} + i \, \Gamma\right)} k_{1}\right)} u_{0}^{2}}{2 \, \Gamma \rho_{0} u_{0}^{2} + {\left(\Gamma^{2} - {\left(\Gamma^{3} - 2 \, \Gamma^{2} + \Gamma\right)} \phi\right)} u_{0}^{3} - 3 \, q_{0}^{2} \rho_{0} - {\left(2 \, \Gamma q_{0}^{2} - \rho_{0}^{2}\right)} u_{0}} & \frac{{\left(3 i \, \Gamma - 5 i\right)} k_{1} q_{0} \rho_{0} u_{0} + {\left(2 i \, \Gamma^{2} - 4 i \, \Gamma\right)} k_{1} q_{0} u_{0}^{2}}{2 \, {\left(2 \, \Gamma \rho_{0} u_{0}^{2} + {\left(\Gamma^{2} - {\left(\Gamma^{3} - 2 \, \Gamma^{2} + \Gamma\right)} \phi\right)} u_{0}^{3} - 3 \, q_{0}^{2} \rho_{0} - {\left(2 \, \Gamma q_{0}^{2} - \rho_{0}^{2}\right)} u_{0}\right)}} & \frac{2 i \, \Gamma k_{1} q_{0} u_{0} + 3 i \, k_{1} q_{0} \rho_{0}}{2 \, {\left(2 \, \Gamma \rho_{0} u_{0}^{2} + {\left(\Gamma^{2} - {\left(\Gamma^{3} - 2 \, \Gamma^{2} + \Gamma\right)} \phi\right)} u_{0}^{3} - 3 \, q_{0}^{2} \rho_{0} - {\left(2 \, \Gamma q_{0}^{2} - \rho_{0}^{2}\right)} u_{0}\right)}} \\ \frac{{\left(i \, \Gamma^{3} - 2 i \, \Gamma^{2} + i \, \Gamma\right)} k_{1} \phi \rho_{0} u_{0}^{3} + {\left(i \, \Gamma^{4} - 2 i \, \Gamma^{3} + i \, \Gamma^{2}\right)} k_{1} \phi u_{0}^{4}}{2 \, \Gamma \rho_{0}^{2} u_{0}^{2} + {\left(\Gamma^{2} - {\left(\Gamma^{3} - 2 \, \Gamma^{2} + \Gamma\right)} \phi\right)} \rho_{0} u_{0}^{3} - 3 \, q_{0}^{2} \rho_{0}^{2} - {\left(2 \, \Gamma q_{0}^{2} \rho_{0} - \rho_{0}^{3}\right)} u_{0}} & -\frac{{\left(i \, \Gamma^{3} - 2 i \, \Gamma^{2} + i \, \Gamma\right)} k_{1} \phi \rho_{0} u_{0}^{2} + {\left(i \, \Gamma^{3} - 2 i \, \Gamma^{2} + i \, \Gamma\right)} k_{1} \phi u_{0}^{3} + {\left(-3 i \, \Gamma + 3 i\right)} k_{1} q_{0}^{2} \rho_{0} + {\left(-2 i \, \Gamma^{2} + 2 i \, \Gamma\right)} k_{1} q_{0}^{2} u_{0}}{2 \, \Gamma \rho_{0} u_{0}^{2} + {\left(\Gamma^{2} - {\left(\Gamma^{3} - 2 \, \Gamma^{2} + \Gamma\right)} \phi\right)} u_{0}^{3} - 3 \, q_{0}^{2} \rho_{0} - {\left(2 \, \Gamma q_{0}^{2} - \rho_{0}^{2}\right)} u_{0}} & \frac{{\left(-5 i \, \Gamma^{2} + i \, \Gamma\right)} k_{1} q_{0} \rho_{0} u_{0}^{2} + 12 i \, k_{1} q_{0}^{3} \rho_{0} + {\left(-2 i \, \Gamma^{3} k_{1} + {\left(4 i \, \Gamma^{3} - 8 i \, \Gamma^{2} + 4 i \, \Gamma\right)} k_{1} \phi\right)} q_{0} u_{0}^{3} + {\left(8 i \, \Gamma k_{1} q_{0}^{3} + {\left(-3 i \, \Gamma + i\right)} k_{1} q_{0} \rho_{0}^{2}\right)} u_{0}}{2 \, {\left(2 \, \Gamma \rho_{0} u_{0}^{2} + {\left(\Gamma^{2} - {\left(\Gamma^{3} - 2 \, \Gamma^{2} + \Gamma\right)} \phi\right)} u_{0}^{3} - 3 \, q_{0}^{2} \rho_{0} - {\left(2 \, \Gamma q_{0}^{2} - \rho_{0}^{2}\right)} u_{0}\right)}} & \frac{-2 i \, \Gamma^{2} k_{1} q_{0} u_{0}^{2} - 5 i \, \Gamma k_{1} q_{0} \rho_{0} u_{0} - 3 i \, k_{1} q_{0} \rho_{0}^{2}}{2 \, {\left(2 \, \Gamma \rho_{0} u_{0}^{2} + {\left(\Gamma^{2} - {\left(\Gamma^{3} - 2 \, \Gamma^{2} + \Gamma\right)} \phi\right)} u_{0}^{3} - 3 \, q_{0}^{2} \rho_{0} - {\left(2 \, \Gamma q_{0}^{2} - \rho_{0}^{2}\right)} u_{0}\right)}} \end{array}\right)$ $\displaystyle \left(\begin{array}{r} \delta_{\rho} \\ \delta_{u} \\ \delta_{u_{1}} \\ \delta_{q} \end{array}\right)$
Analytic eigenvalues and eigenvectors in the $k_1, k_2 \rightarrow 0$ limit :
[($\displaystyle 0$, [$\displaystyle \left(1,\,0,\,0,\,0\right)$, $\displaystyle \left(0,\,1,\,0,\,0\right)$, $\displaystyle \left(0,\,0,\,1,\,0\right)$, $\displaystyle \left(0,\,0,\,0,\,1\right)$], $\displaystyle 4$)]
Numerical eigenvalues and eigenvectors for $k_1 =$ $\displaystyle 1$ , $k_2$ = $\displaystyle 0$ :
-------------------------- ('Eigenvalue = ', 0.2809399731849139*I) (delta_rho, ' = ', 0.9536920617838186) (delta_u, ' = ', 0.13646860180943948 - 6.046083483348376e-17*I) (delta_u1, ' = ', -0.2679302222642111 + 3.357013677842634e-17*I) (delta_q, ' = ', 0.007820997900275337 + 3.011269430191138e-17*I) -------------------------- ('Eigenvalue = ', -1.8778381627448937e-16 + 0.17904252569333254*I) (delta_rho, ' = ', 0.9826577819932445) (delta_u, ' = ', 0.05500422065041492 + 7.28583859910259e-17*I) (delta_u1, ' = ', -0.17593753118027855 - 1.1102230246251565e-16*I) (delta_q, ' = ', 0.020104833273271446 - 5.204170427930421e-18*I) -------------------------- ('Eigenvalue = ', -8.712827791417993e-17 - 0.1725011980833519*I) (delta_rho, ' = ', 0.9844135849570124) (delta_u, ' = ', 0.04081844236148004 - 2.3852447794681098e-17*I) (delta_u1, ' = ', 0.169812522814612 - 1.231653667943533e-16*I) (delta_q, ' = ', -0.020674999651815577 + 2.5370330836160804e-17*I) -------------------------- ('Eigenvalue = ', 1.1516822403160456e-16 - 0.42146374731162467*I) (delta_rho, ' = ', 0.8750616904231109) (delta_u, ' = ', 0.305443974870819 + 2.6242081138666747e-16*I) (delta_u1, ' = ', 0.36880677917456817 + 1.4405483973549372e-16*I) (delta_q, ' = ', 0.07037453945741526 + 1.2741271679096485e-16*I)
eigenvecs_right = M.subs(k1=0, k2=0).eigenvectors_right()

solve(1/eigenvecs_right[0][0] == 0, phi)

[$\displaystyle \phi = \frac{\Gamma^{2} u_{0}^{2} + 2 \, \Gamma \rho_{0} u_{0} + \rho_{0}^{2}}{{\left(\Gamma^{3} - 2 \, \Gamma^{2} + \Gamma\right)} u_{0}^{2}}$]
linearize(gamma/cs**4).subs(delta_u==0, delta_rho==0).simplify_full()

Mk1  = (M - omega*identity_matrix(4))

epsilon = var('epsilon')
dispersion = Mk1.subs(tau=1/epsilon).determinant().numerator().simplify_full()

A = jacobian([Eqn_rho, Eqn_u, Eqn_u1, Eqn_q], [delta_rho_dt, delta_u_dt, delta_u1_dt, delta_q_dt])
B = jacobian([Eqn_rho, Eqn_u, Eqn_u1, Eqn_q], [delta_rho, delta_u, delta_u1, delta_q])
B = B.apply_map(lambda x: x/(I*k1))



M_characteristic = omega*A - B
dispersion_characteristic = M_characteristic.determinant().simplify_full().numerator()
dispersion_characteristic

$\displaystyle 2 \, \Gamma^{3} \omega^{4} \phi u_{0}^{3} + 4 \, \Gamma^{3} \omega^{3} \phi q_{0} u_{0}^{2} - 2 \, \Gamma^{4} \omega^{2} \phi u_{0}^{3} - 4 \, \Gamma^{2} \omega^{4} \phi u_{0}^{3} - 8 \, \Gamma^{2} \omega^{3} \phi q_{0} u_{0}^{2} + 2 \, \Gamma^{3} \omega^{2} \phi \rho_{0} u_{0}^{2} - 2 \, \Gamma^{2} \omega^{4} u_{0}^{3} + 8 \, \Gamma^{3} \omega^{2} \phi u_{0}^{3} + 2 \, \Gamma \omega^{4} \phi u_{0}^{3} + 4 \, \Gamma \omega^{4} q_{0}^{2} u_{0} - 4 \, \Gamma^{2} \omega^{3} q_{0} u_{0}^{2} - 4 \, \Gamma^{3} \omega \phi q_{0} u_{0}^{2} + 4 \, \Gamma \omega^{3} \phi q_{0} u_{0}^{2} - 4 \, \Gamma \omega^{4} \rho_{0} u_{0}^{2} - 4 \, \Gamma^{2} \omega^{2} \phi \rho_{0} u_{0}^{2} + 2 \, \Gamma^{3} \omega^{2} u_{0}^{3} - 2 \, \Gamma^{4} \phi u_{0}^{3} - 10 \, \Gamma^{2} \omega^{2} \phi u_{0}^{3} + 6 \, \omega^{4} q_{0}^{2} \rho_{0} - 2 \, \Gamma \omega^{3} q_{0} \rho_{0} u_{0} - 2 \, \omega^{4} \rho_{0}^{2} u_{0} + 8 \, \Gamma \omega^{3} q_{0} u_{0}^{2} + 8 \, \Gamma^{2} \omega \phi q_{0} u_{0}^{2} + 2 \, \Gamma^{2} \omega^{2} \rho_{0} u_{0}^{2} + 2 \, \Gamma \omega^{2} \phi \rho_{0} u_{0}^{2} - 2 \, \Gamma^{2} \omega^{2} u_{0}^{3} + 6 \, \Gamma^{3} \phi u_{0}^{3} + 4 \, \Gamma \omega^{2} \phi u_{0}^{3} + 3 \, \omega^{3} q_{0} \rho_{0}^{2} - 4 \, \Gamma \omega^{2} q_{0}^{2} u_{0} + 9 \, \omega^{3} q_{0} \rho_{0} u_{0} - 4 \, \Gamma \omega \phi q_{0} u_{0}^{2} - 2 \, \Gamma \omega^{2} \rho_{0} u_{0}^{2} - 6 \, \Gamma^{2} \phi u_{0}^{3} - 6 \, \omega^{2} q_{0}^{2} \rho_{0} - \Gamma \omega q_{0} \rho_{0} u_{0} + 2 \, \Gamma \phi u_{0}^{3} + \omega q_{0} \rho_{0} u_{0}$
epsilon = var('epsilon')
dispersion_characteristic_num = dispersion_characteristic.subs(tau=1/epsilon).subs(epsilon=0, Gamma=4/3, rho0=0, u0=1, phi=12/5, q0=0.088)
dispersion_characteristic_num

$\displaystyle -2.80314311111111 \, \omega^{4} + 0.438044444444444 \, \omega^{3} + 1.61795792592593 \, \omega^{2} - 0.125155555555555 \, \omega - 0.237037037037041$
omega_roots = solve(dispersion_characteristic_num==0, omega, solution_dict=True)
for omega_root in omega_roots:
print omega_root[omega].n(digits=10)

-0.5072419886 - 8.583892882e-14*I -0.4966945818 + 8.280535444e-14*I 0.5503543613 - 1.339117090e-14*I 0.6098512353 + 1.642474527e-14*I
135/16*dispersion_characteristic.subs(tau=1/epsilon).subs(epsilon=0, Gamma=4/3, phi=12/5, rho0=0, u0=1)

$\displaystyle 45 \, \omega^{4} q_{0}^{2} - 24 \, \omega^{4} + 42 \, \omega^{3} q_{0} - 45 \, \omega^{2} q_{0}^{2} + 14 \, \omega^{2} - 12 \, \omega q_{0} - 2$