%auto typeset_mode(True)
balbusaur
A framework for automated linear analysis
Initialize the operators d_dt(), d_dX1() and d_dX2() acting on variables:
Density ρ
Internal energy u
Velocity in X1 direction u1
Velocity in X2 direction u2
Velocity in X3 direction u3
Magnetic field in X1 direction B1
Magnetic field in X2 direction B2
Magnetic field in X3 direction B3
Heat flux magnitude q
Pressure anisotropy magnitude dp
Mean variables : ρ0, u0, u01, u02, u03, B01, B02, B03, q0, dp0
Perturbed variables : δρ, δu, δu1, δu2, δu3, δB1, δB2, δB3, δq, δ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:
EVOLVE_B_FIELDS : 0 or 1
CONDUCTION : 0 or 1
VISCOSITY : 0 or 1
FAKE_EMHD : 0 or 1
TURN_OFF_MEAN_B2 : 0 or 1
TURN_OFF_K2_PERTURBATIONS : 0 or 1
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 = 02Γρ02u02+(Γ2−(Γ3−2Γ2+Γ)ϕ)ρ0u03−3q02ρ02−(2Γq02ρ0−ρ03)u0(2iΓ3−4iΓ2+2iΓ)k1ϕq0u032Γρ02u02+(Γ2−(Γ3−2Γ2+Γ)ϕ)ρ0u03−3q02ρ02−(2Γq02ρ0−ρ03)u0(−iΓ3+2iΓ2−iΓ)k1ϕu032Γρ02u02+(Γ2−(Γ3−2Γ2+Γ)ϕ)ρ0u03−3q02ρ02−(2Γq02ρ0−ρ03)u0(iΓ3−2iΓ2+iΓ)k1ϕρ0u03+(iΓ4−2iΓ3+iΓ2)k1ϕu040−2Γρ0u02+(Γ2−(Γ3−2Γ2+Γ)ϕ)u03−3q02ρ0−(2Γq02−ρ02)u0(−2iΓ+2i)k1q0ρ0u0+((2iΓ3−4iΓ2+2iΓ)k1ϕ+(−2iΓ2+2iΓ)k1)q0u022Γρ0u02+(Γ2−(Γ3−2Γ2+Γ)ϕ)u03−3q02ρ0−(2Γq02−ρ02)u0(−iΓ+i)k1ρ0u0+((iΓ3−2iΓ2+iΓ)k1ϕ+(−iΓ2+iΓ)k1)u02−2Γρ0u02+(Γ2−(Γ3−2Γ2+Γ)ϕ)u03−3q02ρ0−(2Γq02−ρ02)u0(iΓ3−2iΓ2+iΓ)k1ϕρ0u02+(iΓ3−2iΓ2+iΓ)k1ϕu03+(−3iΓ+3i)k1q02ρ0+(−2iΓ2+2iΓ)k1q02u0−ik1ρ02Γρ0u02+(Γ2−(Γ3−2Γ2+Γ)ϕ)u03−3q02ρ0−(2Γq02−ρ02)u0−2iΓ2k1ρ0u03+5ik1q02ρ0u0−(iΓ3k1−(iΓ4−2iΓ3+iΓ2)k1ϕ)u04+(4iΓk1q02−iΓk1ρ02)u022(2Γρ0u02+(Γ2−(Γ3−2Γ2+Γ)ϕ)u03−3q02ρ0−(2Γq02−ρ02)u0)(3iΓ−5i)k1q0ρ0u0+(2iΓ2−4iΓ)k1q0u022(2Γρ0u02+(Γ2−(Γ3−2Γ2+Γ)ϕ)u03−3q02ρ0−(2Γq02−ρ02)u0)(−5iΓ2+iΓ)k1q0ρ0u02+12ik1q03ρ0+(−2iΓ3k1+(4iΓ3−8iΓ2+4iΓ)k1ϕ)q0u03+(8iΓk1q03+(−3iΓ+i)k1q0ρ02)u002Γρ0u02+(Γ2−(Γ3−2Γ2+Γ)ϕ)u03−3q02ρ0−(2Γq02−ρ02)u0−2iΓk1ρ0u02−ik1ρ02u0−(iΓ2k1−(iΓ3−2iΓ2+iΓ)k1ϕ)u032(2Γρ0u02+(Γ2−(Γ3−2Γ2+Γ)ϕ)u03−3q02ρ0−(2Γq02−ρ02)u0)2iΓk1q0u0+3ik1q0ρ02(2Γρ0u02+(Γ2−(Γ3−2Γ2+Γ)ϕ)u03−3q02ρ0−(2Γq02−ρ02)u0)−2iΓ2k1q0u02−5iΓk1q0ρ0u0−3ik1q0ρ02 δρδuδu1δq
Analytic eigenvalues and eigenvectors in the k1,k2→0 limit :
[(0, [(1,0,0,0), (0,1,0,0), (0,0,1,0), (0,0,0,1)], 4)]
Numerical eigenvalues and eigenvectors for k1= 1 , k2 = 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)
[ϕ=(Γ3−2Γ2+Γ)u02Γ2u02+2Γρ0u0+ρ02]
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
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
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