%auto typeset_mode(True)
# Defining the phase space variables var('x, v, t')
(x, v, t)
# Defining the constants var('k, m')
(k, m)
# Definitions of the functions that vary with space and time: var('v_b, T, rho')
(vb, T, ρ)
# Defining collision time-scale var('tau')
τ
# Definition of the distribution function distribution_function_f_0 = rho * (sqrt(m/(2*pi*k*T))) * exp(-m*(v - v_b)^2/(2*k*T)) distribution_function_f_0
21ρπTkme(−2Tkm(v−vb)2)
# Declaring the positive quantities: assume(m>0) assume(k>0) assume(T>0) assume(rho>0)
# Redefining the variables that vary with space and time: v_bulk = function('v_{bulk}')(x, t) T = function('T')(x, t) rho = function('rho')(x, t)
# Defining the various velocity moments: zeroth_moment = integrate(distribution_function_f_0, v, -infinity, infinity).subs(T=T, rho=rho, v_b=v_bulk) first_moment = integrate(distribution_function_f_0*v, v, -infinity, infinity).subs(T=T, rho=rho, v_b=v_bulk) second_moment = integrate(distribution_function_f_0*v**2, v, -infinity, infinity).subs(T=T, rho=rho, v_b=v_bulk) third_moment = integrate(distribution_function_f_0*v**3, v, -infinity, infinity).subs(T=T, rho=rho, v_b=v_bulk) fourth_moment = integrate(distribution_function_f_0*v**4, v, -infinity, infinity).subs(T=T, rho=rho, v_b=v_bulk)
# Simplifying the expressions zeroth_moment = zeroth_moment.expand().simplify() first_moment = first_moment.expand() second_moment = second_moment.expand() third_moment = third_moment.expand() fourth_moment = fourth_moment.expand()
# Obtaining the Euler equations euler_continuity_eqn = zeroth_moment.diff(t) + first_moment.diff(x) euler_momentum_eqn = first_moment.diff(t) + second_moment.diff(x) euler_energy_eqn = second_moment.diff(t) + third_moment.diff(x)
# Obtaining the NS equations # The expression for the distribution function f_approx is f = f_0 - tau * v * df_0/dx # Due to mass, momentum and energy conservation: # \int f dv = \int f_0 dv # \int v f dv = \int v f_0 dv # \int v^2 f dv = \int v^2 f_0 dv NS_continuity_eqn = zeroth_moment.diff(t) + first_moment.diff(x) NS_momentum_eqn = first_moment.diff(t) + second_moment.diff(x) NS_energy_eqn = second_moment.diff(t) + (third_moment - tau * fourth_moment.diff(x)).diff(x)
# Linearizing the system: var('rho0, delta_rho, T0, delta_T, v0, delta_v, I, k1, delta_rho_dt, delta_T_dt, delta_v_dt')
(ρ0, δρ, T0, δT, v0, δv, I, k1, δρdt, δTdt, δvdt)
def linearize(term): return taylor(term, (delta_rho, 0), \ (delta_v, 0), \ (delta_T, 0), \ (delta_rho_dt, 0), \ (delta_v_dt, 0), \ (delta_T_dt, 0), 1 \ ).simplify_full() eqn_continuity = linearize(NS_continuity_eqn.subs(rho==rho0+delta_rho, rho.diff(x)==I*k1*delta_rho, rho.diff(t)==delta_rho_dt, \ T==T0+delta_T, T.diff(x)==I*k1*delta_T, T.diff(t)==delta_T_dt, \ v_bulk==v0+delta_v, v_bulk.diff(x)==I*k1*delta_v, v_bulk.diff(t)==delta_v_dt ) ) eqn_momentum = linearize(NS_momentum_eqn.subs(rho==rho0+delta_rho, rho.diff(x)==I*k1*delta_rho, rho.diff(t)==delta_rho_dt, \ T==T0+delta_T, T.diff(x)==I*k1*delta_T, T.diff(t)==delta_T_dt, \ v_bulk==v0+delta_v, v_bulk.diff(x)==I*k1*delta_v, v_bulk.diff(t)==delta_v_dt ) ) eqn_energy = linearize(NS_energy_eqn.subs(rho==rho0+delta_rho, rho.diff(x)==I*k1*delta_rho, rho.diff(t)==delta_rho_dt, rho.diff(x, 2)==-k1**2*delta_rho,\ T==T0+delta_T, T.diff(x)==I*k1*delta_T, T.diff(t)==delta_T_dt, T.diff(x, 2)==-k1**2*delta_T,\ v_bulk==v0+delta_v, v_bulk.diff(x)==I*k1*delta_v, v_bulk.diff(t)==delta_v_dt, v_bulk.diff(x, 2)==-k1**2*delta_v ) ) eqns = [eqn_continuity == 0, eqn_momentum == 0, eqn_energy == 0] delta_vars = [delta_rho, delta_v, delta_T] delta_vars_dt = [delta_rho_dt, delta_v_dt, delta_T_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())
Linearized system :
δρdtδvdtδTdt = −Ik1v0−mρ0IT0kk1−kmρ0k12m2τv04+6T0kk12mτv02+3T02k2k12τ−Ik1ρ0−Ik1v0−k2(2k12mτv03+6T0kk12τv0+IT0kk1)0−mIkk1−m6k12mτv02+6T0kk12τ+Ik1mv0 δρδvδT
rho0_num = 1 v0_num = 0 T0_num = 1 k1_num = 2 * pi tau_num = 0.01 M_numerical = M.subs(rho0=rho0_num, v0=v0_num, T0=T0_num, k1=k1_num, tau = tau_num) M_numerical = M_numerical.change_ring(CDF) eigenvecs = M_numerical.eigenvectors_right()