48 lines
1.0 KiB
Python
48 lines
1.0 KiB
Python
|
#! /usr/bin/env python3
|
||
|
|
import sympy as sp
|
||
|
|
|
||
|
|
sp.init_printing(use_unicode=True, use_latex=True)
|
||
|
|
|
||
|
|
rho = sp.symbols("rho", real=True, positive=True, nonzero=True)
|
||
|
|
ru, rv, e = sp.symbols("rhou rhov e", real=True)
|
||
|
|
gamma = sp.symbols("gamma", real=True, positive=True, nonzero=True)
|
||
|
|
p = (gamma - 1) * (e - (ru ** 2 + rv ** 2) / (2 * rho))
|
||
|
|
c = sp.sqrt(gamma * p / rho)
|
||
|
|
|
||
|
|
u = ru / rho
|
||
|
|
v = rv / rho
|
||
|
|
|
||
|
|
E = sp.Matrix([rho * u, rho * u * u + p, rho * u * v, u * (e + p)])
|
||
|
|
F = sp.Matrix([rho * v, rho * u * v, rho * v * v + p, v * (e + p)])
|
||
|
|
|
||
|
|
A = E.jacobian([rho, ru, rv, e])
|
||
|
|
B = F.jacobian([rho, ru, rv, e])
|
||
|
|
|
||
|
|
print()
|
||
|
|
for key in A.eigenvals():
|
||
|
|
sp.pprint(key)
|
||
|
|
print()
|
||
|
|
for key in B.eigenvals():
|
||
|
|
sp.pprint(key)
|
||
|
|
|
||
|
|
# Eigenvalues are u, u +- c and v, v +- c
|
||
|
|
# sp.pprint((key - c).simplify()) # Remainder after subtracting c
|
||
|
|
|
||
|
|
# kx, ky = sp.symbols("kx ky")
|
||
|
|
# Arot = kx*A + ky*B
|
||
|
|
|
||
|
|
S, L = A.diagonalize()
|
||
|
|
|
||
|
|
sp.pprint(S)
|
||
|
|
sp.pprint(L)
|
||
|
|
|
||
|
|
Labs = (L * L) ** 0.5
|
||
|
|
|
||
|
|
plus = S * (L + Labs) * S.inv() / 2
|
||
|
|
sp.pprint(plus)
|
||
|
|
minus = S * (L - Labs) * S.inv() / 2
|
||
|
|
sp.pprint(minus)
|
||
|
|
|
||
|
|
S = plus - minus
|
||
|
|
sp.pprint(S)
|