This module is currently under development
This module is organized as follows
We will now use the cmdline_helper.py NRPy+ module to create the source directory within the IllinoisGRMHD
NRPy+ directory, if it does not exist yet.
# Step 0: Creation of the IllinoisGRMHD source directory
# Step 0a: Add NRPy's directory to the path
# https://stackoverflow.com/questions/16780014/import-file-from-parent-directory
import os,sys
nrpy_dir_path = os.path.join("..","..")
if nrpy_dir_path not in sys.path:
sys.path.append(nrpy_dir_path)
# Step 0b: Load up cmdline_helper and create the directory
import cmdline_helper as cmd
IGM_src_dir_path = os.path.join("..","src")
cmd.mkdir(IGM_src_dir_path)
# Step 0c: Create the output file path
outfile_path__eigen__C = os.path.join(IGM_src_dir_path,"eigen.C")
In this tutorial notebook we will implement an algorithm to evaluate the eigenvalues of a $3\times3$ symmetric matrix. Our method will be analytical and will follow closely this discussion.
Let $\mathcal{M}$ be a $3\times3$ symmetric matrix,
$$ \mathcal{M} = \begin{pmatrix} M_{11} & M_{12} & M_{13}\\ M_{12} & M_{22} & M_{23}\\ M_{13} & M_{23} & M_{33} \end{pmatrix}\ . $$To obtain the eigenvalues of $\mathcal{M}$, we must solve the characteristic equation
$$ \det\left(\lambda I_{3\times3} - \mathcal{M}\right) = 0\ , $$where $\lambda$ represents the eigenvalues of $\mathcal{M}$ and $I_{3\times3} = {\rm diag}\left(1,1,1\right)$ is the $3\times3$ identity matrix. For this particular case, the characteristic equation of $\mathcal{M}$ is then given by
$$ \lambda^{3} - {\rm tr}\left(\mathcal{M}\right)\lambda^{2} + \left[\frac{{\rm tr}\left(\mathcal{M}^{2}\right) - {\rm tr}\left(\mathcal{M}\right)^{2}}{2}\right]\lambda - \det\left(\mathcal{M}\right) = 0\ . $$Now let $\mathcal{M} = n\mathcal{N} + mI_{3\times3}$, so that the matrices $\mathcal{M}$ and $\mathcal{N}$ have the same eigenvectors. Then, $\kappa$ is an eigenvalue of $\mathcal{N}$ if, and only if, $\lambda = n\kappa + m$ is an eigenvalue of $\mathcal{M}$. Now, let us look at the following identities:
$$ \mathcal{N} = \frac{1}{n}\left(\mathcal{M} - mI_{3\times3}\right)\ . $$Choosing $m \equiv \frac{1}{3}{\rm tr}\left(\mathcal{M}\right)$, we get
$$ {\rm tr}\left(\mathcal{N}\right) = \frac{1}{n}\left(\mathcal{M} - 3m\right)=0\ . $$Also,
$$ {\rm tr}\left(\mathcal{N}^{2}\right) = \frac{1}{n^{2}}\left[N_{11}^{2}+N_{22}^{2}+N_{33}^{2}+2\left(N_{12}^{2}+N_{13}^{2}+N_{23}^{2}\right)\right]\ , $$so that if we choose $n\equiv\sqrt{\frac{N_{11}^{2}+N_{22}^{2}+N_{33}^{2}+2\left(N_{12}^{2}+N_{13}^{2}+N_{23}^{2}\right)}{6}}$ we get
$$ {\rm tr}\left(\mathcal{N}^{2}\right) = 6\ . $$Then, if we look at the characteristic equation for the matrix $\mathcal{N}$,
$$ \kappa^{3} - {\rm tr}\left(\mathcal{N}\right)\kappa^{2} + \left[\frac{{\rm tr}\left(\mathcal{N}^{2}\right) - {\rm tr}\left(\mathcal{N}\right)^{2}}{2}\right]\kappa - \det\left(\mathcal{N}\right) = 0\ , $$we see that it can be greatly simplified with our choices of $m$ and $n$,
$$ \kappa^{3} - 3\kappa - \det\left(\mathcal{N}\right) = 0\ . $$Further simplification of this characteristic equation can be obtained by using
$$ \begin{align} \kappa &\equiv 2\cos\phi\ ,\\ \cos\left(3\phi\right) &= 4\cos^{3}\phi - 3\cos\phi\ , \end{align} $$so that
$$ \begin{align} 0 &= 8\cos^{3}\phi - 6\cos\phi - \det\left(\mathcal{N}\right)\\ &= 2\cos\left(3\phi\right) - \det\left(\mathcal{N}\right)\\ \implies \phi &= \frac{1}{3}\arccos\frac{\det\left(\mathcal{N}\right)}{2} + \frac{2k\pi}{3}\ ,\ k=0,1,2\ , \end{align} $$which, finally, yields
$$ \boxed{\kappa\left(k\right) = 2\cos\left(\frac{1}{3}\arccos\frac{\det\left(\mathcal{N}\right)}{2}+\frac{2k\pi}{3}\right)}\ . $$Once we have $\kappa$, we can find the eigenvectors of $\mathcal{M}$ doing
$$ \boxed{ \begin{align} \lambda_{1} &= m + 2n\kappa(0)\\ \lambda_{2} &= m + 2n\kappa(1)\\ \lambda_{3} &= 3m - \lambda_{1} - \lambda_{2} \end{align} }\ , $$where we have used the fact that ${\rm tr}\left(\mathcal{M}\right)=\lambda_{1}+\lambda_{2}+\lambda_{3}$ to compute $\lambda_{3}$.
eigen.C
[Back to top]¶eigen.C
[Back to top]¶In the algorithm below, we define the following quantities
$$ \boxed{ \begin{align} \mathcal{K} &= \mathcal{M} - mI_{3\times3}\\ m &= \frac{{\rm tr}\left(\mathcal{M}\right)}{3}\\ q &= \frac{\det\left(\mathcal{K}\right)}{2}\\ p &= n^{2} = \frac{{\rm tr}\left(\mathcal{K}^{2}\right)}{6} \end{align} }\ . $$We these definitions, we have the following quantities to be implemented:
$$ \boxed{ m = \frac{\left(M_{11} + M_{22} + M_{33}\right)}{3} }\ . $$The matrix $\mathcal{K}$ is simply
$$ \boxed{ \mathcal{K} = \begin{pmatrix} M_{11}-m & M_{12} & M_{13}\\ M_{12} & M_{22}-m & M_{23}\\ M_{13} & M_{23} & M_{33}-m \end{pmatrix} }\ . $$Straightforwardly, we have
$$ \boxed{q = \frac{K_{11}K_{22}K_{33} + K_{12}K_{23}K_{13} + K_{13}K_{12}K_{23} - K_{13}K_{22}K_{13} - K_{12}K_{12}K_{33} - K_{11}K_{23}K_{23} }{2} }\ . $$Since $\mathcal{K}$ is symmetric as well, we have
$$ \boxed{p = \frac{K_{11}^{2} + K_{22}^{2} + K_{33}^{2} + 2\left(K_{12}^{2} + K_{13}^{2} + K_{23}^{2}\right)}{6}}\ . $$%%writefile $outfile_path__eigen__C
//
// This subroutine calcualtes the eigenvalues of a real, symmetric 3x3
// matrix M={{M11,M12,M13},{M12,M22,M23},{M13,M23,M33}} based on the
// algorithm described in
// http://en.wikipedia.org/wiki/Eigenvalue_algorithm#Eigenvalues_of_3.C3.973_matrices
// which simply solve the cubic equation Det( M - lamnda I)=0 analytically.
// The eigenvalues are stored in lam1, lam2 and lam3.
//
void eigenvalues_3by3_real_sym_matrix(CCTK_REAL & lam1, CCTK_REAL & lam2, CCTK_REAL & lam3,
CCTK_REAL M11, CCTK_REAL M12, CCTK_REAL M13, CCTK_REAL M22, CCTK_REAL M23, CCTK_REAL M33)
{
CCTK_REAL m = (M11 + M22 + M33)/3.0;
CCTK_REAL K11 = M11 - m, K12 = M12, K13 = M13, K22 = M22-m, K23 = M23, K33=M33-m;
CCTK_REAL q = 0.5* (K11*K22*K33 + K12*K23*K13 + K13*K12*K23 - K13*K22*K13
- K12*K12*K33 - K11*K23*K23);
CCTK_REAL p = ( SQR(K11) + SQR(K22) + SQR(K33) + 2.0*(SQR(K12) + SQR(K13) + SQR(K23) ) )/6.0;
Writing ../src/eigen.C
%%writefile -a $outfile_path__eigen__C
CCTK_REAL phi;
CCTK_REAL p32 = sqrt(p*p*p);
if (fabs(q) >= fabs(p32) ) {
phi = 0.0;
} else {
phi = acos(q/p32)/3.0;
}
if (phi<0.0) phi += M_PI/3.0;
Appending to ../src/eigen.C
Finally, the eigenvalues are computed using
$$ \boxed{ \begin{align} \lambda_{1} &= m + 2\sqrt{p}\cos\phi\\ \lambda_{2} &= m - \sqrt{p}\cos\phi - \sqrt{3p}\sin\phi\\ \lambda_{3} &= m - \sqrt{p}\cos\phi + \sqrt{3p}\sin\phi \end{align} }\ . $$%%writefile -a $outfile_path__eigen__C
CCTK_REAL sqrtp = sqrt(p);
CCTK_REAL sqrtp_cosphi = sqrtp*cos(phi);
CCTK_REAL sqrtp_sqrt3_sinphi = sqrtp*sqrt(3.0)*sin(phi);
lam1 = m + 2.0*sqrtp_cosphi;
lam2 = m - sqrtp_cosphi - sqrtp_sqrt3_sinphi;
lam3 = m - sqrtp_cosphi + sqrtp_sqrt3_sinphi;
}
Appending to ../src/eigen.C
# Verify if the code generated by this tutorial module
# matches the original IllinoisGRMHD source code
# First download the original IllinoisGRMHD source code
import urllib
from os import path
original_IGM_file_url = "https://bitbucket.org/zach_etienne/wvuthorns/raw/5611b2f0b17135538c9d9d17c7da062abe0401b6/IllinoisGRMHD/src/eigen.C"
original_IGM_file_name = "eigen-original.C"
original_IGM_file_path = os.path.join(IGM_src_dir_path,original_IGM_file_name)
# Then download the original IllinoisGRMHD source code
# We try it here in a couple of ways in an attempt to keep
# the code more portable
try:
original_IGM_file_code = urllib.request.urlopen(original_IGM_file_url).read().decode("utf-8")
# Write down the file the original IllinoisGRMHD source code
with open(original_IGM_file_path,"w") as file:
file.write(original_IGM_file_code)
except:
try:
original_IGM_file_code = urllib.urlopen(original_IGM_file_url).read().decode("utf-8")
# Write down the file the original IllinoisGRMHD source code
with open(original_IGM_file_path,"w") as file:
file.write(original_IGM_file_code)
except:
# If all else fails, hope wget does the job
!wget -O $original_IGM_file_path $original_IGM_file_url
# Perform validation
Validation__eigen__C = !diff $original_IGM_file_path $outfile_path__eigen__C
if Validation__eigen__C == []:
# If the validation passes, we do not need to store the original IGM source code file
!rm $original_IGM_file_path
print("Validation test for eigen.C: PASSED!")
else:
# If the validation fails, we keep the original IGM source code file
print("Validation test for eigen.C: FAILED!")
# We also print out the difference between the code generated
# in this tutorial module and the original IGM source code
print("Diff:")
for diff_line in Validation__eigen__C:
print(diff_line)
Validation test for eigen.C: FAILED! Diff: 16a17 > 24a26 > 31a34 >
The following code cell converts this Jupyter notebook into a proper, clickable $\LaTeX$-formatted PDF file. After the cell is successfully run, the generated PDF may be found in the root NRPy+ tutorial directory, with filename Tutorial-IllinoisGRMHD__eigen.pdf (Note that clicking on this link may not work; you may need to open the PDF file through another means).
latex_nrpy_style_path = os.path.join(nrpy_dir_path,"latex_nrpy_style.tplx")
#!jupyter nbconvert --to latex --template $latex_nrpy_style_path --log-level='WARN' Tutorial-IllinoisGRMHD__eigen.ipynb
#!pdflatex -interaction=batchmode Tutorial-IllinoisGRMHD__eigen.tex
#!pdflatex -interaction=batchmode Tutorial-IllinoisGRMHD__eigen.tex
#!pdflatex -interaction=batchmode Tutorial-IllinoisGRMHD__eigen.tex
!rm -f Tut*.out Tut*.aux Tut*.log