CoCalc Public Filesgrim / src / problem / linear_modes / Balbusaur.sagewsOpen with one click!
Authors: Mani Chandra, Francois Foucart
Views : 47
Compute Environment: Ubuntu 18.04 (Deprecated)
%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) = (00ik1ρ00(2iΓ34iΓ2+2iΓ)k1ϕq0u032Γρ02u02+(Γ2(Γ32Γ2+Γ)ϕ)ρ0u033q02ρ02(2Γq02ρ0ρ03)u0(2iΓ+2i)k1q0ρ0u0+((2iΓ34iΓ2+2iΓ)k1ϕ+(2iΓ2+2iΓ)k1)q0u022Γρ0u02+(Γ2(Γ32Γ2+Γ)ϕ)u033q02ρ0(2Γq02ρ02)u02iΓ2k1ρ0u03+5ik1q02ρ0u0(iΓ3k1(iΓ42iΓ3+iΓ2)k1ϕ)u04+(4iΓk1q02iΓk1ρ02)u022Γρ0u02+(Γ2(Γ32Γ2+Γ)ϕ)u033q02ρ0(2Γq02ρ02)u02iΓk1ρ0u02ik1ρ02u0(iΓ2k1(iΓ32iΓ2+iΓ)k1ϕ)u032Γρ0u02+(Γ2(Γ32Γ2+Γ)ϕ)u033q02ρ0(2Γq02ρ02)u0(iΓ3+2iΓ2iΓ)k1ϕu032Γρ02u02+(Γ2(Γ32Γ2+Γ)ϕ)ρ0u033q02ρ02(2Γq02ρ0ρ03)u0(iΓ+i)k1ρ0u0+((iΓ32iΓ2+iΓ)k1ϕ+(iΓ2+iΓ)k1)u022Γρ0u02+(Γ2(Γ32Γ2+Γ)ϕ)u033q02ρ0(2Γq02ρ02)u0(3iΓ5i)k1q0ρ0u0+(2iΓ24iΓ)k1q0u022(2Γρ0u02+(Γ2(Γ32Γ2+Γ)ϕ)u033q02ρ0(2Γq02ρ02)u0)2iΓk1q0u0+3ik1q0ρ02(2Γρ0u02+(Γ2(Γ32Γ2+Γ)ϕ)u033q02ρ0(2Γq02ρ02)u0)(iΓ32iΓ2+iΓ)k1ϕρ0u03+(iΓ42iΓ3+iΓ2)k1ϕu042Γρ02u02+(Γ2(Γ32Γ2+Γ)ϕ)ρ0u033q02ρ02(2Γq02ρ0ρ03)u0(iΓ32iΓ2+iΓ)k1ϕρ0u02+(iΓ32iΓ2+iΓ)k1ϕu03+(3iΓ+3i)k1q02ρ0+(2iΓ2+2iΓ)k1q02u02Γρ0u02+(Γ2(Γ32Γ2+Γ)ϕ)u033q02ρ0(2Γq02ρ02)u0(5iΓ2+iΓ)k1q0ρ0u02+12ik1q03ρ0+(2iΓ3k1+(4iΓ38iΓ2+4iΓ)k1ϕ)q0u03+(8iΓk1q03+(3iΓ+i)k1q0ρ02)u02(2Γρ0u02+(Γ2(Γ32Γ2+Γ)ϕ)u033q02ρ0(2Γq02ρ02)u0)2iΓ2k1q0u025iΓk1q0ρ0u03ik1q0ρ022(2Γρ0u02+(Γ2(Γ32Γ2+Γ)ϕ)u033q02ρ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,k20k_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(Γ32Γ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ϕq0u022Γ4ω2ϕu034Γ2ω4ϕu038Γ2ω3ϕq0u02+2Γ3ω2ϕρ0u022Γ2ω4u03+8Γ3ω2ϕu03+2Γω4ϕu03+4Γω4q02u04Γ2ω3q0u024Γ3ωϕq0u02+4Γω3ϕq0u024Γω4ρ0u024Γ2ω2ϕρ0u02+2Γ3ω2u032Γ4ϕu0310Γ2ω2ϕu03+6ω4q02ρ02Γω3q0ρ0u02ω4ρ02u0+8Γω3q0u02+8Γ2ωϕq0u02+2Γ2ω2ρ0u02+2Γω2ϕρ0u022Γ2ω2u03+6Γ3ϕu03+4Γω2ϕu03+3ω3q0ρ024Γω2q02u0+9ω3q0ρ0u04Γωϕq0u022Γω2ρ0u026Γ2ϕu036ω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ω20.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ω4q0224ω4+42ω3q045ω2q02+14ω212ωq02\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