| Hosted by CoCalc | Download
%auto
typeset_mode(True)

balbusaur\mathtt{balbusaur}

A framework for automated linear analysis

Initialize the operators d_dt()\mathtt{d\_dt()}, d_dX1()\mathtt{d\_dX1()} and d_dX2()\mathtt{d\_dX2()} acting on variables:

  1. Density ρ\rho

  2. Internal energy uu

  3. Velocity in X1X^1 direction u1u^1

  4. Velocity in X2X^2 direction u2u^2

  5. Velocity in X3X^3 direction u3u^3

  6. Magnetic field in X1X^1 direction B1B^1

  7. Magnetic field in X2X^2 direction B2B^2

  8. Magnetic field in X3X^3 direction B3B^3

  9. Heat flux magnitude qq

  10. Pressure anisotropy magnitude dpdp

  • Mean variables : ρ0\rho_0, u0u_0, u01u^1_0, u02u^2_0, u03u^3_0, B01B^1_0, B02B^2_0, B03B^3_0, q0q_0, dp0dp_0

  • Perturbed variables : δρ\delta_\rho, δu\delta_u, δu1\delta_{u^1}, δu2\delta_{u^2}, δu3\delta_{u^3}, δB1\delta_{B^1}, δB2\delta_{B^2}, δB3\delta_{B^3}, δq\delta_q, δdp\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. EVOLVE_B_FIELDS\mathtt{EVOLVE\_B\_FIELDS} : 0 or 1

  2. CONDUCTION\mathtt{CONDUCTION} : 0 or 1

  3. VISCOSITY\mathtt{VISCOSITY} : 0 or 1

  4. FAKE_EMHD\mathtt{FAKE\_EMHD} : 0 or 1

  5. TURN_OFF_MEAN_B2\mathtt{TURN\_OFF\_MEAN\_B2} : 0 or 1

  6. TURN_OFF_K2_PERTURBATIONS\mathtt{TURN\_OFF\_K2\_PERTURBATIONS} : 0 or 1

  7. PRINCIPAL_COEFFICIENTS\mathtt{PRINCIPAL\_COEFFICIENTS} : 0 or 1

# Spatiotemporal variables
t, omega, k1, k2 = var('t, omega, k1, k2')

# Constants:
# Gamma : Adiabatic index
# 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 :
(δρdtδudtδu1dtδqdt)\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) = (00−i k1ρ00(2i Γ3−4i Γ2+2i Γ)k1ϕq0u032 Γρ02u02+(Γ2−(Γ3−2 Γ2+Γ)ϕ)ρ0u03−3 q02ρ02−(2 Γq02ρ0−ρ03)u0−(−2i Γ+2i)k1q0ρ0u0+((2i Γ3−4i Γ2+2i Γ)k1ϕ+(−2i Γ2+2i Γ)k1)q0u022 Γρ0u02+(Γ2−(Γ3−2 Γ2+Γ)ϕ)u03−3 q02ρ0−(2 Γq02−ρ02)u0−2i Γ2k1ρ0u03+5i k1q02ρ0u0−(i Γ3k1−(i Γ4−2i Γ3+i Γ2)k1ϕ)u04+(4i Γk1q02−i Γk1ρ02)u022 Γρ0u02+(Γ2−(Γ3−2 Γ2+Γ)ϕ)u03−3 q02ρ0−(2 Γq02−ρ02)u0−2i Γk1ρ0u02−i k1ρ02u0−(i Γ2k1−(i Γ3−2i Γ2+i Γ)k1ϕ)u032 Γρ0u02+(Γ2−(Γ3−2 Γ2+Γ)ϕ)u03−3 q02ρ0−(2 Γq02−ρ02)u0(−i Γ3+2i Γ2−i Γ)k1ϕu032 Γρ02u02+(Γ2−(Γ3−2 Γ2+Γ)ϕ)ρ0u03−3 q02ρ02−(2 Γq02ρ0−ρ03)u0(−i Γ+i)k1ρ0u0+((i Γ3−2i Γ2+i Γ)k1ϕ+(−i Γ2+i Γ)k1)u022 Γρ0u02+(Γ2−(Γ3−2 Γ2+Γ)ϕ)u03−3 q02ρ0−(2 Γq02−ρ02)u0(3i Γ−5i)k1q0ρ0u0+(2i Γ2−4i Γ)k1q0u022 (2 Γρ0u02+(Γ2−(Γ3−2 Γ2+Γ)ϕ)u03−3 q02ρ0−(2 Γq02−ρ02)u0)2i Γk1q0u0+3i k1q0ρ02 (2 Γρ0u02+(Γ2−(Γ3−2 Γ2+Γ)ϕ)u03−3 q02ρ0−(2 Γq02−ρ02)u0)(i Γ3−2i Γ2+i Γ)k1ϕρ0u03+(i Γ4−2i Γ3+i Γ2)k1ϕu042 Γρ02u02+(Γ2−(Γ3−2 Γ2+Γ)ϕ)ρ0u03−3 q02ρ02−(2 Γq02ρ0−ρ03)u0−(i Γ3−2i Γ2+i Γ)k1ϕρ0u02+(i Γ3−2i Γ2+i Γ)k1ϕu03+(−3i Γ+3i)k1q02ρ0+(−2i Γ2+2i Γ)k1q02u02 Γρ0u02+(Γ2−(Γ3−2 Γ2+Γ)ϕ)u03−3 q02ρ0−(2 Γq02−ρ02)u0(−5i Γ2+i Γ)k1q0ρ0u02+12i k1q03ρ0+(−2i Γ3k1+(4i Γ3−8i Γ2+4i Γ)k1ϕ)q0u03+(8i Γk1q03+(−3i Γ+i)k1q0ρ02)u02 (2 Γρ0u02+(Γ2−(Γ3−2 Γ2+Γ)ϕ)u03−3 q02ρ0−(2 Γq02−ρ02)u0)−2i Γ2k1q0u02−5i Γk1q0ρ0u0−3i k1q0ρ022 (2 Γρ0u02+(Γ2−(Γ3−2 Γ2+Γ)ϕ)u03−3 q02ρ0−(2 Γq02−ρ02)u0))\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) (δρδuδu1δq)\displaystyle \left(\begin{array}{r} \delta_{\rho} \\ \delta_{u} \\ \delta_{u_{1}} \\ \delta_{q} \end{array}\right)
Analytic eigenvalues and eigenvectors in the k1,k2→0k_1, k_2 \rightarrow 0 limit :
[(0\displaystyle 0, [(1, 0, 0, 0)\displaystyle \left(1,\,0,\,0,\,0\right), (0, 1, 0, 0)\displaystyle \left(0,\,1,\,0,\,0\right), (0, 0, 1, 0)\displaystyle \left(0,\,0,\,1,\,0\right), (0, 0, 0, 1)\displaystyle \left(0,\,0,\,0,\,1\right)], 4\displaystyle 4)]
Numerical eigenvalues and eigenvectors for k1=k_1 = 1\displaystyle 1 , k2k_2 = 0\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)
[ϕ=Γ2u02+2 Γρ0u0+ρ02(Γ3−2 Γ2+Γ)u02\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
2 Γ3ω4ϕu03+4 Γ3ω3ϕq0u02−2 Γ4ω2ϕu03−4 Γ2ω4ϕu03−8 Γ2ω3ϕq0u02+2 Γ3ω2ϕρ0u02−2 Γ2ω4u03+8 Γ3ω2ϕu03+2 Γω4ϕu03+4 Γω4q02u0−4 Γ2ω3q0u02−4 Γ3ωϕq0u02+4 Γω3ϕq0u02−4 Γω4ρ0u02−4 Γ2ω2ϕρ0u02+2 Γ3ω2u03−2 Γ4ϕu03−10 Γ2ω2ϕu03+6 ω4q02ρ0−2 Γω3q0ρ0u0−2 ω4ρ02u0+8 Γω3q0u02+8 Γ2ωϕq0u02+2 Γ2ω2ρ0u02+2 Γω2ϕρ0u02−2 Γ2ω2u03+6 Γ3ϕu03+4 Γω2ϕu03+3 ω3q0ρ02−4 Γω2q02u0+9 ω3q0ρ0u0−4 Γωϕq0u02−2 Γω2ρ0u02−6 Γ2ϕu03−6 ω2q02ρ0−Γωq0ρ0u0+2 Γϕu03+ωq0ρ0u0\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
−2.80314311111111 ω4+0.438044444444444 ω3+1.61795792592593 ω2−0.125155555555555 ω−0.237037037037041\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)
45 ω4q02−24 ω4+42 ω3q0−45 ω2q02+14 ω2−12 ωq0−2\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