

In [1]:

    
from sympy import init_session
init_session()









    



IPython console for SymPy 1.0 (Python 2.7.12-64-bit) (ground types: gmpy)

These commands were executed:
>>> from __future__ import division
>>> 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 http://docs.sympy.org/1.0/

Euler Equations

The Euler equations in primitive variable form, $q = (\rho, u, p)^\intercal$ appear as:

$$q_t + A(q) q_x = 0$$

with the matrix $A(q)$:

$$A(q) = \left ( \begin{array}{ccc} u & \rho & 0 \\ 0 & u & 1/\rho \\ 0 & \gamma p & u \end{array} \right ) $$

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.



In [2]:

    
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









    Out[2]:





$$\left[\begin{matrix}u & \rho & 0\\0 & u & \frac{1}{\rho}\\0 & c^{2} \rho & u\end{matrix}\right]$$

The eigenvalues are the speeds at which information propagates with. SymPy returns them as a dictionary, giving the multiplicity for each eigenvalue.



In [3]:

    
A.eigenvals()









    Out[3]:





$$\left \{ u : 1, \quad - c + u : 1, \quad c + u : 1\right \}$$

The right eigenvectors are what SymPy gives natively. For a given eigenvalue, $\lambda$, these satisfy:

$$A r = \lambda r$$

Right Eigenvectors



In [4]:

    
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



In [5]:

    
r[0]









    Out[5]:





$$\left[\begin{matrix}1\\0\\0\end{matrix}\right]$$

this corresponds to the eigenvalue



In [6]:

    
lam[0]









    Out[6]:





$$u$$

1-st right eigenvector



In [7]:

    
r[1]









    Out[7]:





$$\left[\begin{matrix}1\\- \frac{c}{\rho}\\c^{2}\end{matrix}\right]$$

this corresponds to the eigenvalue



In [8]:

    
lam[1]









    Out[8]:





$$- c + u$$

2-nd right eigenvector



In [9]:

    
r[2]









    Out[9]:





$$\left[\begin{matrix}1\\\frac{c}{\rho}\\c^{2}\end{matrix}\right]$$

this corresponds to the eigenvalue



In [10]:

    
lam[2]









    Out[10]:





$$c + u$$

Here they are as a matrix, $R$, in order from smallest to largest eigenvalue



In [11]:

    
R = zeros(3,3)
R[:,0] = r[1]
R[:,1] = r[0]
R[:,2] = r[2]
R









    Out[11]:





$$\left[\begin{matrix}1 & 1 & 1\\- \frac{c}{\rho} & 0 & \frac{c}{\rho}\\c^{2} & 0 & c^{2}\end{matrix}\right]$$

Left Eigenvectors

The left eigenvectors satisfy:

$$l A = \lambda l$$

SymPy doesn't have a method to get left eigenvectors directly, so we take the transpose of this expression:

$$(l A)^\intercal = A^\intercal l^\intercal = \lambda l^\intercal$$

Therefore, the transpose of the left eigenvectors, $l^\intercal$, are the right eigenvectors of transpose of $A$



In [12]:

    
B = A.transpose()
B









    Out[12]:





$$\left[\begin{matrix}u & 0 & 0\\\rho & u & c^{2} \rho\\0 & \frac{1}{\rho} & u\end{matrix}\right]$$



In [13]:

    
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}$



In [14]:

    
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



In [15]:

    
l[0]









    Out[15]:





$$\left[\begin{matrix}1 & 0 & - \frac{1}{c^{2}}\end{matrix}\right]$$

1-st left eigenvector



In [16]:

    
l[1]









    Out[16]:





$$\left[\begin{matrix}0 & - \frac{\rho}{2 c} & \frac{1}{2 c^{2}}\end{matrix}\right]$$

2-nd left eigenvector



In [17]:

    
l[2]









    Out[17]:





$$\left[\begin{matrix}0 & \frac{\rho}{2 c} & \frac{1}{2 c^{2}}\end{matrix}\right]$$



In [ ]:

Entropy formulation

here we write the system in terms of $q_s = (\rho, u, s)^\intercal$, where the system is

$${q_s}_t + A_s(q_s) {q_s}_x = 0$$

and

$$ A_s = \left (\begin{matrix}u & \rho & 0\\ \frac{c^{2}}{\rho} & u & \frac{p_{s}}{\rho}\\ 0 & 0 & u\end{matrix}\right) $$



In [18]:

    
ps = symbols('p_s')

As = Matrix([[u, rho, 0], [c**2/rho, u, ps/rho], [0, 0, u]])
As









    Out[18]:





$$\left[\begin{matrix}u & \rho & 0\\\frac{c^{2}}{\rho} & u & \frac{p_{s}}{\rho}\\0 & 0 & u\end{matrix}\right]$$



In [19]:

    
As.eigenvals()









    Out[19]:





$$\left \{ u : 1, \quad - c + u : 1, \quad c + u : 1\right \}$$



In [20]:

    
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]



In [21]:

    
r[0], lam[0]









    Out[21]:





$$\left ( \left[\begin{matrix}1\\0\\- \frac{c^{2}}{p_{s}}\end{matrix}\right], \quad u\right )$$



In [22]:

    
r[1], lam[1]









    Out[22]:





$$\left ( \left[\begin{matrix}1\\- \frac{c}{\rho}\\0\end{matrix}\right], \quad - c + u\right )$$



In [23]:

    
r[2], lam[2]









    Out[23]:





$$\left ( \left[\begin{matrix}1\\\frac{c}{\rho}\\0\end{matrix}\right], \quad c + u\right )$$

left eigenvectors



In [24]:

    
Bs = As.transpose()
L = B.eigenvects()
l = []
laml = []
for (ev, _, ltmp) in L:
    l.append(ltmp[0].transpose())
    laml.append(ev)

normalization



In [25]:

    
for n in range(len(l)):
    if lam[n] == laml[n]:
        ltmp = l[n]
        p = ltmp.dot(r[n])
        l[n] = ltmp/p



In [26]:

    
simplify(l[0])









    Out[26]:





$$\left[\begin{matrix}\frac{p_{s}}{p_{s} + 1} & 0 & - \frac{p_{s}}{c^{2} \left(p_{s} + 1\right)}\end{matrix}\right]$$



In [27]:

    
l[1]









    Out[27]:





$$\left[\begin{matrix}0 & - \frac{\rho}{c} & \frac{1}{c^{2}}\end{matrix}\right]$$



In [28]:

    
l[2]









    Out[28]:





$$\left[\begin{matrix}0 & \frac{\rho}{c} & \frac{1}{c^{2}}\end{matrix}\right]$$

2-d system



In [29]:

    
rho, u, v, c = symbols('rho u v c')

A = Matrix([[u, rho, 0, 0], [0, u, 0, rho**-1], [0,0, u, 0], [0, c**2 * rho, 0, u]])
A









    Out[29]:





$$\left[\begin{matrix}u & \rho & 0 & 0\\0 & u & 0 & \frac{1}{\rho}\\0 & 0 & u & 0\\0 & c^{2} \rho & 0 & u\end{matrix}\right]$$



In [30]:

    
A.eigenvals()









    Out[30]:





$$\left \{ u : 2, \quad - c + u : 1, \quad c + u : 1\right \}$$



In [31]:

    
R = A.eigenvects()   # this returns a tuple for each eigenvector with multiplicity -- unpack it
r = []
lam = []
for (ev, _, rtmp) in R:
    for rv in rtmp:
        r.append(rv)
        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]
    if not v[0] == 0:
        r[n] = v/v[0]



In [32]:

    
r[0], lam[0]









    Out[32]:





$$\left ( \left[\begin{matrix}1\\0\\0\\0\end{matrix}\right], \quad u\right )$$



In [33]:

    
r[1], lam[1]









    Out[33]:





$$\left ( \left[\begin{matrix}0\\0\\1\\0\end{matrix}\right], \quad u\right )$$



In [34]:

    
r[2], lam[2]









    Out[34]:





$$\left ( \left[\begin{matrix}1\\- \frac{c}{\rho}\\0\\c^{2}\end{matrix}\right], \quad - c + u\right )$$



In [35]:

    
r[3], lam[3]









    Out[35]:





$$\left ( \left[\begin{matrix}1\\\frac{c}{\rho}\\0\\c^{2}\end{matrix}\right], \quad c + u\right )$$



In [ ]: