from sympy import init_session
init_session(use_latex="mathjax")
IPython console for SymPy 1.13.3 (Python 3.10.15-64-bit) (ground types: python)
These commands were executed:
>>> from sympy import *
>>> x, y, z, t = symbols('x y z t')
>>> k, m, n = symbols('k m n', integer=True)
>>> f, g, h = symbols('f g h', cls=Function)
>>> init_printing()
Documentation can be found at https://docs.sympy.org/1.13.3/
Euler equations eigensystem#
Primitive variable form#
The Euler equations in primitive variable form, \({\bf q} = (\rho, u, p)^\intercal\) appear as:
with the matrix \({\bf A}({\bf q})\):
The sound speed is related to the adiabatic index, \(\gamma\), as \(c^2 = \gamma p /\rho\).
We can represent this matrix symbolically in SymPy and explore its eigensystem.
from sympy.abc import rho
rho, u, c = symbols('rho u c')
A = Matrix([[u, rho, 0], [0, u, rho**-1], [0, c**2 * rho, u]])
A
The eigenvalues are the speeds at which information propagates with. SymPy returns them as a dictionary, giving the multiplicity for each eigenvalue.
A.eigenvals()
The right eigenvectors are what SymPy gives natively. For a given eigenvalue, \(\lambda\), these satisfy:
Right Eigenvectors#
R = A.eigenvects() # this returns a tuple for each eigenvector with multiplicity -- unpack it
r = []
lam = []
for (ev, _, rtmp) in R:
r.append(rtmp[0])
lam.append(ev)
# we can normalize them anyway we want, so let's make the first entry 1
for n in range(len(r)):
v = r[n]
r[n] = v/v[0]
0-th right eigenvector#
r[0]
this corresponds to the eigenvalue
lam[0]
1-st right eigenvector#
r[1]
this corresponds to the eigenvalue
lam[1]
2-nd right eigenvector#
r[2]
this corresponds to the eigenvalue
lam[2]
Here they are as a matrix, \({\bf R}\), in order from smallest to largest eigenvalue
R = zeros(3,3)
R[:,0] = r[1]
R[:,1] = r[0]
R[:,2] = r[2]
R
Left Eigenvectors#
The left eigenvectors satisfy:
We’ll find these by taking the transpose of \({\bf A}\)
Therefore, the transpose of the left eigenvectors, \({\bf l}^\intercal\), are the right eigenvectors of transpose of \({\bf A}\)
B = A.transpose()
B
L = B.eigenvects()
l = []
laml = []
for (ev, _, ltmp) in L:
l.append(ltmp[0].transpose())
laml.append(ev)
Traditionally, we normalize these such that \(l^{(\mu)} \cdot r^{(\nu)} = \delta_{\mu\nu}\)
for n in range(len(l)):
if lam[n] == laml[n]:
ltmp = l[n]
p = ltmp.dot(r[n])
l[n] = ltmp/p
0-th left eigenvector#
l[0]
1-st left eigenvector#
l[1]
2-nd left eigenvector#
l[2]
Entropy formulation#
here we write the system in terms of \({\bf q}_s = (\rho, u, s)^\intercal\), where the system is
and
ps = symbols('p_s')
As = Matrix([[u, rho, 0], [c**2/rho, u, ps/rho], [0, 0, u]])
As
As.eigenvals()
R = As.eigenvects() # this returns a tuple for each eigenvector with multiplicity -- unpack it
r = []
lam = []
for (ev, _, rtmp) in R:
r.append(rtmp[0])
lam.append(ev)
# we can normalize them anyway we want, so let's make the first entry 1
for n in range(len(r)):
v = r[n]
r[n] = v/v[0]
r[0]
this corresponds to eigenvalue
lam[0]
r[1]
this corresponds to eigenvalue
lam[1]
r[2]
this corresponds to eigenvalue
lam[2]
left eigenvectors#
Bs = As.transpose()
L = B.eigenvects()
l = []
laml = []
for (ev, _, ltmp) in L:
l.append(ltmp[0].transpose())
laml.append(ev)
normalization
for n in range(len(l)):
if lam[n] == laml[n]:
ltmp = l[n]
p = ltmp.dot(r[n])
l[n] = ltmp/p
simplify(l[0])
l[1]
l[2]