(%i1) /* working with Jameson, "Eigenvalues, Eigenvectors, and Symmetrization of the Magneto-Hydrodynamic (MHD) Equations" 1996

- shifting the pressure conservative variable to the last
- sorting the eigenvalues from least to greatest
- verifying that the derivative of conservative wrt symmetrized variables and inverse are in fact inverses of one another
- verifying that the eigenvectors of symmetrized matrix are orthonormal
- verifying that the reconstructed system is the flux

NOTICE I was using eigenvectors() until I found out about this ...
http://stackoverflow.com/questions/19135430/diagonalize-a-symmetric-matrix-in-maxima
... then I started using ev(similaritytransform(M), hermitianmatrix=true) instead
*/"";

Result

(%i2) load("diag");

Result

(%i3) load("invert");

Result

(%i4) load("eigen");

Result

(%i5) swapsides(eqn) := rhs(eqn) = lhs(eqn)$

(%i6) makeEigenMatrices(eigenvectorResults) :=
block([n, reps, evs],
    n : length(eigenvectorResults[1][1]),
    evs : eigenvectorResults[1][1],
    reps : eigenvectorResults[1][2],
    [diag(apply(append, makelist(makelist(evs[j],i,1,reps[j]),j,1,n))),
        transpose(apply(matrix, apply(append, eigenvectorResults[2])))])$

(%i7) /* normalize columns of a matrix */
normalizeMatrix(m) := m . diag(makelist(1/sqrt(col(m,i).col(m,i)),i,1,length(m)))$

(%i8) assume(mu>0, rho>0, c>=0, ca>=0, cs>=0, cf>=0);

Result

(%i9) xs : [x,y,z];

Result

(%i10) vs : [vx,vy,vz];

Result

(%i11) Bs : [Bx,By,Bz];

Result

(%i12) /* conservative variables */
us : [u1,u2,u3,u4,u5,u6,u7,u8];

Result

(%i13) /* primitive variables */
ws : [rho,vx,vy,vz,Bx,By,Bz,p];

Result

(%i14) depends(append(ws,[P,Z,E,c,ca,cs,cf]),append(xs,[t]));

Result

(%i15) depends(us, append(xs,[t]));

Result

(%i16) vSq : vs.vs;

Result

(%i17) BSq : Bs.Bs;

Result

(%i18) divB : sum('diff(Bs[i],xs[i]),i,1,3);

Result

(%i19) VdotB : sum(Bs[i]*vs[i],i,1,3);

Result

(%i20) /* Z = total specific energy (hydro + magnetic)
Z = E + B^2 / (2 rho mu)
Z = e + 1/2 v^2 + B^2 / (2 rho mu) */
Z_from_E_B_rho : Z = E + BSq / (2 * rho * mu);

Result

(%i21) /* total pressure (hydro + magnetic) */
P_from_p_B : P = p + BSq / (2 * mu);

Result

(%i22) p_from_P_B : solve(P_from_p_B, p)[1];

Result

(%i23) /* total specific hydro energy from internal specific energy + kinetic specific energy */
E_from_e_v : E = e + vSq/2;

Result

(%i24) e_from_E_v : -(E_from_e_v - e - E);

Result

(%i25) /* fluid pressure for perfect gas.
p = fluid pressure.
e = internal specific hydro energy
 */
p_from_rho_e : p = (%gamma - 1)*rho*e;

Result

(%i26) p_from_rho_e / ((%gamma - 1) * rho)$
rhs(%) = lhs(%)$
e_from_p_rho : %;

Result

(%i29) /* p = fluid pressure for perfect gas
   E = total specific energy (hydro)
   v^2/2 = kinetic specific energy
   E - v^2/2 = e = internal specific energy
*/
p_from_rho_E_v : subst([e_from_E_v], p_from_rho_e);

Result

(%i30) solve(p_from_rho_E_v, E)[1]$
E_from_rho_v_p : %;

Result

(%i32) subst([E_from_rho_v_p], Z_from_E_B_rho)$
ratsimp(%)$
Z_from_rho_v_p_B : %;

Result

(%i35) /* speed of sound for perfect gas */
cSq_from_rho_p : c^2 = %gamma * p / rho;

Result

(%i36) p_from_rho_c : solve(cSq_from_rho_p, p)[1];

Result

(%i37) /* internal specific enthalpy */
h_from_e : h = %gamma * e;

Result

(%i38) h_from_rho_p : subst([e_from_p_rho], h_from_e);

Result

(%i39) p_from_rho_h : solve(h_from_rho_p, p)[1];

Result

(%i40) /* total specific enthlapy */
H_from_h_v : H = h + vSq/2;

Result

(%i41) H_from_rho_v_p : ratsimp(subst([h_from_rho_p], H_from_h_v));

Result

(%i42) H_from_v_c : ratsimp(subst([p_from_rho_c],H_from_rho_v_p));

Result

(%i43) subst([p_from_rho_c], H_from_rho_v_p)$
solve(%, c^2)[1]$
cSq_from_v_H : %;

Result

(%i46) subst([H_from_h_v],cSq_from_v_H)$
ratsimp(%)$
cSq_from_h : %;

Result

(%i49) caSq_from_rho_B : ca^2 = Bx^2 / (rho * mu);

Result

(%i50) /* is the Alfven speed based on |Bx| or Bx?
Maxima says Bx, but it is paired up with Bx and -Bx, which is some permutation of |Bx| and -|Bx|*/
ca_from_rho_B : sqrt(caSq_from_rho_B);

Result

(%i51) cStar^2 = 1/2*(c^2 + BSq/(rho*mu))$
ratsimp(%)$
cStarSq_from_rho_B_c : %;

Result

(%i54) cStarSq_from_rho_B_c^2$
cStarFourth_from_rho_B_c : %;

Result

(%i56) cf^2 = cStar^2 + sqrt(cStar^4 - c^2 * ca^2)$
cfSq_from_c_ca_cStar : %;

Result

(%i58) subst([cStarSq_from_rho_B_c, cStarFourth_from_rho_B_c, caSq_from_rho_B], cfSq_from_c_ca_cStar)$
ratsimp(%)$
cfSq_from_rho_B_c : %;

Result

(%i61) sqrt(cfSq_from_rho_B_c)$
ratsimp(%)$
cf_from_rho_B_c : %;

Result

(%i64) cs^2 = cStar^2 - sqrt(cStar^4 - c^2 * ca^2)$
csSq_from_c_ca_cStar : %;

Result

(%i66) subst([cStarSq_from_rho_B_c, cStarFourth_from_rho_B_c, caSq_from_rho_B], csSq_from_c_ca_cStar)$
ratsimp(%)$
csSq_from_rho_B_c : %;

Result

(%i69) sqrt(csSq_from_rho_B_c)$
ratsimp(%)$
cs_from_rho_B_c : %;

Result

(%i72) /* relation between conservative and primitive */
[u1=rho, u2=rho*vx, u3=rho*vy, u4=rho*vz, u5=Bx, u6=By, u7=Bz, u8=rho*Z]$
subst([Z_from_rho_v_p_B], %)$
ratsimp(%)$
u_for_w : %;

Result

(%i76) solve(u_for_w, ws)[1]$
ratsimp(%)$
w_for_u : %;

Result

(%i79) /* du/dw */
makelist(rhs(u_for_w[i]),i,1,8)$
jacobian(%, ws)$
ratsimp(%)$
dCons_dPrim : %;

Result

(%i83) invert(dCons_dPrim)$
ratsimp(%)$
dPrim_dCons : %;

Result

(%i86) /* verify inverse */
is(ratsimp(dCons_dPrim.dPrim_dCons)=ident(8));
is(ratsimp(dPrim_dCons.dCons_dPrim)=ident(8));

Result

(%i88) /* equations of motion, used for flux */"";

Result

(%i89) continuity_eqn : 'diff(rho,t) + sum('diff(rho*vs[i],xs[i]),i,1,3) = 0;

Result

(%i90) makelist(
    'diff(rho*vs[i],t)
    + sum(
        'diff(rho*vs[i]*vs[j]
        - Bs[i]*Bs[j]/mu
        + kron_delta(i,j)*P,xs[j])
    ,j,1,3)
    = -Bs[i]/mu*divB
,i,1,3);
subst([P_from_p_B], %)$
ratsimp(%)$
momentum_eqns : %;

Result

(%i94) magnetic_field_eqns : makelist(
    'diff(Bs[i],t)
    + sum(
        'diff(
            vs[j]*Bs[i]-vs[i]*Bs[j]
        ,xs[j])
    ,j,1,3)
    = -vs[i]*divB
,i,1,3);

Result

(%i95) 'diff(rho*Z,t) + sum('diff((rho*Z+P)*vs[j]-VdotB*Bs[j]/mu,xs[j]),j,1,3) = -VdotB / mu * divB$
subst([Z_from_rho_v_p_B], %)$
subst([P_from_p_B],%) /* the Jameson paper says p, but every single other paper says P */$
ratsimp(%)$
energy_eqn : %$

(%i100) append([continuity_eqn], momentum_eqns, magnetic_field_eqns, [energy_eqn])$
all_eqns_from_prims : %$

(%i102) subst([w_for_u], all_eqns_from_prims)$
ratsimp(%)$
all_eqns_from_cons : %$

(%i105) all_eqns_from_cons$
expand(ev(%, diff))$
subst([makelist('diff(us[i],y)=0,i,1,8)],%) /* remove y and z derivatives while we're here ... we're only looking at the x-axis */ $
subst([makelist('diff(us[i],z)=0,i,1,8)],%)$
all_eqns_from_cons_expanded : %$

(%i110) /* Why can Maxima defer differentiation but can't defer matrix multiplies?
Until then, here's me picking the matrices out of the equations, one by one.

flux along x axis, coeffs wrt conservative variables, neglecting source term (lhs of the eqn excluding time derivative)
du/dt + dF/du du/dx = S du/dx
this is the dF/du
*/
genmatrix(lambda([i,j], coeff(lhs(all_eqns_from_cons_expanded[i]), 'diff(us[j],x))),8,8)$
ratsimp(%)$
fluxDerivWrtCons_from_cons : %$

(%i113) fluxDerivWrtCons_from_cons$
subst([u_for_w], %)$
ratsimp(%)$
fluxDerivWrtCons_from_prims : %$

(%i117) /* Source term along x axis, coeffs wrt conservative variables (rhs of equation).
du/dt + dF/du du/dx = S du/dx
this is the S
*/
genmatrix(lambda([i,j], coeff(rhs(all_eqns_from_cons_expanded[i]), 'diff(us[j],x))),8,8)$
subst([u_for_w], %)$
ratsimp(%)$
sourceDerivWrtCons_from_prims : %;

Result

(%i121) /* du/dt + dF/du du/dx = dS/du du/dx
du/dt + (dF/du - S) du/dx = 0
A = dF/du - S */
fluxDerivWrtCons_from_prims - sourceDerivWrtCons_from_prims$
subst([u_for_w], %)$
ratsimp(%)$
quasiLinearDerivWrtCons_for_prims : %;

Result

(%i125) /*
du/dt + A du/dx = 0 ... change in variables
du/dw dw/dt + A du/dw dw/dx = 0 ... lhs multiply by [du/dw]^-1 = dw/du
dw/dt + (dw/du A du/dw) dw/dx = 0
let ATilde = dw/du A du/dw
*/
dPrim_dCons . quasiLinearDerivWrtCons_for_prims . dCons_dPrim$
subst([p_from_rho_c], %)$
ratsimp(%)$
quasiLinearDerivWrtPrims : %;

Result

(%i129) /* verification: A = du/dw ATilde dw/du */
dCons_dPrim . quasiLinearDerivWrtPrims . dPrim_dCons - quasiLinearDerivWrtCons_for_prims$
subst([cSq_from_rho_p],%)$
ratsimp(%)$
is(%=0*ident(8));

Result

(%i133) /* now for a final change in variables to symmetrize A
define the transformation rather than the coordinate chart (could be anholonomic but I'm too lazy ot check) */
dSym_dPrim : matrix(
[0,0,0,0,0,0,0,1/(rho*c)],
[0,1,0,0,0,0,0,0],
[0,0,1,0,0,0,0,0],
[0,0,0,1,0,0,0,0],
[0,0,0,0,1/sqrt(mu*rho),0,0,0],
[0,0,0,0,0,1/sqrt(mu*rho),0,0],
[0,0,0,0,0,0,1/sqrt(mu*rho),0],
[-c/rho,0,0,0,0,0,0,1/(rho*c)]
) * rho/c;

Result

(%i134) dPrim_dSym : invert(dSym_dPrim);

Result

(%i135) /* verify inverse */
is(dPrim_dSym.dSym_dPrim = ident(8));
is(dSym_dPrim.dPrim_dSym = ident(8));

Result

(%i137) /*
du/dt + A du/dx = 0 ... change in variables
du/dw dw/dt + A du/dw dw/dx = 0 ... lhs multiply by [du/dw]^-1 = dw/du
dw/dt + (dw/du A du/dw) dw/dx = 0

... change of variables from primitives (w) to symmetrizing basis (y)
dw/dy dy/dt + (dw/du A du/dw) dw/dy dy/dx = 0 ... right multiply by [dw/dy]^-1 = dy/dw
dy/dt + (dy/dw dw/du A du/dw dw/dy) dy/dx = 0

let ABar = dy/dw ATilde dw/dy = dy/dw dw/du A du/dw dw/dy = dy/du A du/dy
*/
dSym_dPrim . quasiLinearDerivWrtPrims . dPrim_dSym$
ratsimp(%)$
quasiLinearDerivWrtSym : %;

Result

(%i140) is(quasiLinearDerivWrtSym = transpose(quasiLinearDerivWrtSym));

Result

(%i141) /* Note:
let ABar = dy/dw ATilde dw/dy = dy/dw dw/du A du/dw dw/dy = dy/du A du/dy
therefore A = [dy/du]^-1 ABar [du/dy]^-1 = du/dy ABar dy/du = du/dw dw/dy ABar dy/dw dw/du.
So after finding the eigensystem of ABar, transform it by du/dw.dw/dy
*/"";

Result

(%i142) /* du/dy = dCons/dSym = dCons/dPrim . dPrim/dSym */
dCons_dPrim . dPrim_dSym$
subst([cSq_from_v_H],%)$
ratsimp(%)$
dCons_dSym : %;

Result

(%i146) /* Maxima coercion */
block([A,B],
    A : (%gamma-1)/c^2*(vSq-H),
    B : subst([H_from_rho_v_p, p_from_rho_c], A),
    B : ratsimp(B),
    B = A)$
tmpEquality_dSym_dCons_8_1 : %;

Result

(%i148) /* dy/du = dSym/dCons = dSym/dPrim . dPrim/dCons */
block([A],
    A : dSym_dPrim . dPrim_dCons,
    A : ratsimp(A),
    A : subst([tmpEquality_dSym_dCons_8_1],A))$
dSym_dCons : %;

Result

(%i150) /* verification of the dCons/dSym basis */
dCons_dSym.dSym_dCons$
subst([H_from_rho_v_p,p_from_rho_c],%)$
ratsimp(%)$
is(%=ident(8));

Result

(%i154) /* verify we can restore the original flux */
dCons_dSym . quasiLinearDerivWrtSym . dSym_dCons - quasiLinearDerivWrtCons_for_prims$
ratsimp(%)$
subst([H_from_rho_v_p,p_from_rho_c],%)$
ratsimp(%)$
is(%=0*ident(8));

Result

(%i159) /* now for the eigensystem decomposition */
quasiLinearDerivWrtSym_eigenvector_results : ev(similaritytransform(quasiLinearDerivWrtSym), hermitianmatrix=true)$

(%i160) [quasiLinearDerivWrtSym_maximaEigenvalues_notReduced, quasiLinearDerivWrtSym_maximaEigenvectors_notReduced]
: makeEigenMatrices(quasiLinearDerivWrtSym_eigenvector_results)$

(%i161) /* clean up the eigenvalues matrix using the slow, fast, and Alfvens wave speeds */"";

Result

(%i162) cStarSq_from_rho_B_c$
% * 4 * rho * mu$
ratsimp(%)$
rhs(%)=lhs(%)$
tmp1 : %;

Result

(%i167) (cStarFourth_from_rho_B_c - cStar^4) * 4 * mu^2 * rho^2$
solve(%, c^4)[1]$
tmp2 : %;

Result

(%i170) caSq_from_rho_B$
solve(%, Bx^2)[1]$
tmp3 : %;

Result

(%i173) csSq_from_c_ca_cStar$
rhs(%)=lhs(%)$
tmp4 : %;

Result

(%i176) cfSq_from_c_ca_cStar$
rhs(%)=lhs(%)$
tmp5 : %;

Result

(%i179) quasiLinearDerivWrtSym_maximaEigenvalues_notReduced$
subst([tmp1, tmp2, tmp3], %)$
factor(%)$
subst([tmp4, tmp5], %)$
expand(%)$
quasiLinearDerivWrtSym_maximaEigenvalues_reduced : %;

Result

(%i185) /* also, permuting eigenvectors: V_mine = P V_original P^-1 <=> V_original = P^-1 V_mine P
... A = U_original V_original U_original^-1 = U_original P^-1 P V_original P^-1 P U_original^-1
... A = (U_original P^-1) V_mine (U_original P^-1)^-1 = U_mine V_mine U_mine^-1
so U_mine = U_original P^-1
*/"";

Result

(%i186) permuteEigenvalues : matrix(
[0,0,1,0,0,0,0,0],
[0,0,0,0,1,0,0,0],
[1,0,0,0,0,0,0,0],
[0,0,0,0,0,0,1,0],
[0,0,0,0,0,0,0,1],
[0,1,0,0,0,0,0,0],
[0,0,0,0,0,1,0,0],
[0,0,0,1,0,0,0,0]
);

Result

(%i187) permuteEigenvalues . quasiLinearDerivWrtSym_maximaEigenvalues_reduced . invert(permuteEigenvalues)$
quasiLinearDerivWrtSym_eigenvalues_reduced : %;

Result

(%i189) /* clean up the eigenvector matrix using slow, fast, and Alfven wave speeds */"";

Result

(%i190) cStarSq_from_rho_B_c$
% * 2 * rho * mu$
% - Bx^2 - By^2 - Bz^2$
ratsimp(%)$
rhs(%) = lhs(%)$
tmp1b : %;

Result

(%i196) cfSq_from_c_ca_cStar - cStar^2$
rhs(%) = lhs(%)$
tmp6 : %;

Result

(%i199) cStar^2 - csSq_from_c_ca_cStar$
rhs(%) = lhs(%)$
tmp7 : %;

Result

(%i202) is(lhs(tmp6)=lhs(tmp7));
rhs(tmp7)=rhs(tmp6)$
(% + cStar^2 + cs^2) / 2$
tmp8 : %;

Result

(%i206) tmp8$
% * 2$
subst([cStarSq_from_rho_B_c], %)$
% * rho * mu$
expand(%)$
tmp9 : %;

Result

(%i212) tmp9 - BSq - cf^2*rho*mu$
rhs(%) = lhs(%)$
tmp10 : %;

Result

(%i215) tmp9 - BSq - cs^2*rho*mu$
rhs(%) = lhs(%)$
tmp11 : %;

Result

(%i218) quasiLinearDerivWrtSym_maximaEigenvectors_notReduced . invert(permuteEigenvalues)$
quasiLinearDerivWrtSym_eigenvectors_notReduced : %$

(%i220) block([m],
    m : quasiLinearDerivWrtSym_eigenvectors_notReduced,
    for i thru 8 do
        m : ratsimp(m . diag(makelist(if (m[i,j]=0) then 1 else denom(m[i,j]),j,1,8))),
    m)$
subst([tmp1b, tmp2, tmp3],%)$
factor(%)$
subst([tmp4, tmp5, tmp6],%)$
factor(%)$
% . diag([sqrt(mu*rho/2),1/abs(Bz),sqrt(mu*rho/2),1,1,sqrt(mu*rho/2),1/abs(Bz),sqrt(mu*rho/2)])$
subst([tmp1b, caSq_from_rho_B],%)$
% . diag([1/(2*mu^2*rho^2),1,1/(2*mu^2*rho^2),1,1,1/(2*mu^2*rho^2),1,1/(2*mu^2*rho^2)])$
expand(%)$
subst([swapsides(cfSq_from_c_ca_cStar + csSq_from_c_ca_cStar)/2],%)$
expand(%)$
quasiLinearDerivWrtSym_eigenvectors_reduced :%;

Result

(%i232) /* eigenvector orthogonality test. make sure it's A^T A since we're right-multiplying A by a scale matrix */
transpose(quasiLinearDerivWrtSym_eigenvectors_reduced) . quasiLinearDerivWrtSym_eigenvectors_reduced$
expand(%)$
subst([cfSq_from_c_ca_cStar, cf*cfSq_from_c_ca_cStar, cfSq_from_c_ca_cStar^2, cfSq_from_c_ca_cStar^3], %)$
subst([csSq_from_c_ca_cStar, cs*csSq_from_c_ca_cStar, csSq_from_c_ca_cStar^2, csSq_from_c_ca_cStar^3], %)$
subst([caSq_from_rho_B, caSq_from_rho_B^2], %)$
ratsimp(%)$
expand(%)$
subst([cStarSq_from_rho_B_c],%)$
expand(%)$
/* subst([cStarFourth_from_rho_B_c],%) */
%;

Result

(%i242) /* eigenvectors of degenerate cases for Bx=0; By=Bz=0; and Bx=By=Bz=0 */"";

Result

(%i243) /* DEGENERACY BX=BY=BZ=0 */"";

Result

(%i244) quasiLinearDerivWrtSym$
subst([Bx=0,By=0,Bz=0],%)$
quasiLinearDerivWrtSym_zeroBxyz : %;

Result

(%i247) quasiLinearDerivWrtSym_zeroBxyz_eigenvector_results : ev(similaritytransform(quasiLinearDerivWrtSym_zeroBxyz), hermitianmatrix=true);

Result

(%i248) block([A,B],
    [A,B] : makeEigenMatrices(quasiLinearDerivWrtSym_zeroBxyz_eigenvector_results),
    B : B . diag([sqrt(2),sqrt(2),1,1,1,1,1,1]),
    [quasiLinearDerivWrtSym_zeroBxyz_maximaEigenvalues_reduced, quasiLinearDerivWrtSym_zeroBxyz_maximaEigenvectors_reduced] : [expand(A),B]);

Result

(%i249) /* eigenvector orthogonality */
quasiLinearDerivWrtSym_zeroBxyz_maximaEigenvectors_reduced$
transpose(%) . %$
normalizeMatrix(%)$
is(% = ident(8));

Result

(%i253) permuteEigenvalues_zeroBxyz : matrix(
[1,0,0,0,0,0,0,0],
[0,0,1,0,0,0,0,0],
[0,0,0,1,0,0,0,0],
[0,0,0,0,1,0,0,0],
[0,0,0,0,0,0,0,1],
[0,0,0,0,0,1,0,0],
[0,0,0,0,0,0,1,0],
[0,1,0,0,0,0,0,0]);

Result

(%i254) permuteEigenvalues_zeroBxyz . quasiLinearDerivWrtSym_zeroBxyz_maximaEigenvalues_reduced . invert(permuteEigenvalues_zeroBxyz)$
quasiLinearDerivWrtSym_zeroBxyz_eigenvalues_reduced : %;

Result

(%i256) quasiLinearDerivWrtSym_zeroBxyz_maximaEigenvectors_reduced . invert(permuteEigenvalues_zeroBxyz)$
quasiLinearDerivWrtSym_zeroBxyz_eigenvectors_reduced : %;

Result

(%i258) /* eigenvector orthogonality */
quasiLinearDerivWrtSym_zeroBxyz_eigenvectors_reduced$
transpose(%) . %$
% : normalizeMatrix(%)$
is(% = ident(8));

Result

(%i262) /* DEGENERACY BX=0
therefore ca=0, cs=0 */"";

Result

(%i263) ca_from_rho_B$
subst([Bx=0],%)$
ca_zeroBx:%;

Result

(%i266) cStarSq_from_rho_B_c$
subst([Bx=0],%)$
cStarSq_from_rho_B_c_zeroBx : %;

Result

(%i269) cfSq_from_c_ca_cStar$
subst([ca_zeroBx, cStarSq_from_rho_B_c_zeroBx], %)$
% - (By^2 + Bz^2) / (rho * mu)$
ratsimp(%)$
swapsides(%)$
cSq_from_rho_B_cf_zeroBx : %;

Result

(%i275) cSq_from_rho_B_cf_zeroBx$
% + (By^2 + Bz^2)/(rho*mu)$
expand(%)$
swapsides(%)$
cfSq_from_rho_B_c_zeroBx : %;

Result

(%i280) quasiLinearDerivWrtSym$
subst([Bx=0],%)$
quasiLinearDerivWrtSym_zeroBx : %;

Result

(%i283) ev(similaritytransform(quasiLinearDerivWrtSym_zeroBx), hermitianmatrix=true)$
quasiLinearDerivWrtSym_zeroBx_eigenvector_results : %;

Result

(%i285) block([A,B],
    [A,B] : makeEigenMatrices(quasiLinearDerivWrtSym_zeroBx_eigenvector_results),
    A : subst([cSq_from_rho_B_cf_zeroBx], expand(A)),
    B : subst([cSq_from_rho_B_cf_zeroBx], B),
    [quasiLinearDerivWrtSym_zeroBx_maximaEigenvalues_reduced, quasiLinearDerivWrtSym_zeroBx_maximaEigenvectors_reduced] : [A, B]
);

Result

(%i286) /* TODO eigenvector orthogonality */
quasiLinearDerivWrtSym_zeroBx_maximaEigenvectors_reduced$
transpose(%) . %
/*normalizeMatrix(%)
is(% = ident(8))*/$

(%i288) permuteEigenvalues_zeroBx : matrix(
[1,0,0,0,0,0,0,0],
[0,0,0,1,0,0,0,0],
[0,0,0,0,1,0,0,0],
[0,0,0,0,0,1,0,0],
[0,0,0,0,0,0,0,1],
[0,0,0,0,0,0,1,0],
[0,0,1,0,0,0,0,0],
[0,1,0,0,0,0,0,0]);

Result

(%i289) is(permuteEigenvalues_zeroBx . invert(permuteEigenvalues_zeroBx) = ident(8));

Result

(%i290) permuteEigenvalues_zeroBx . quasiLinearDerivWrtSym_zeroBx_maximaEigenvalues_reduced . invert(permuteEigenvalues_zeroBx)$
quasiLinearDerivWrtSym_zeroBx_eigenvalues_reduced : %;

Result

(%i292) quasiLinearDerivWrtSym_zeroBx_maximaEigenvectors_reduced . invert(permuteEigenvalues_zeroBx)$
% . diag([
    sqrt(2*(cf^2*rho*mu-By^2-Bz^2)+2*By^2+2*Bz^2)/sqrt(rho*mu),
    1,
    1,
    1,
    1,
    cf*sqrt(cf^2*rho*mu-By^2)/sqrt(rho*mu),
    -c*Bz*sqrt(cf^2*rho*mu-By^2)/(c*abs(Bz)*sqrt(rho*mu)),
    sqrt(2*(cf^2*rho*mu-By^2-Bz^2)+2*By^2+2*Bz^2)/sqrt(rho*mu)])$
quasiLinearDerivWrtSym_zeroBx_eigenvectors_reduced : %;

Result

(%i295) /* notice that the Bx=0 case matches up with the Bx=By=Bz=0 case except
for the 5th column of the Bx=0 case -- which is composed of [Bz,-By] where the Bx=By=Bz=0 case has [1,0]
there's nothing else in that column, so an extra test will have to be performed to check tangent magnitude
and substitute [1,0] in place of none. But only the 5th row -- for the sake of orthogonality. The 1st and 8th row [By,Bz] should remain zero.*/"";

Result

(%i296) /* eigenvector orthogonality */
transpose(quasiLinearDerivWrtSym_zeroBx_eigenvectors_reduced) . quasiLinearDerivWrtSym_zeroBx_eigenvectors_reduced$
ratsimp(%)$
factor(%)$
subst([cfSq_from_rho_B_c_zeroBx, cfSq_from_rho_B_c_zeroBx^2],%)$
ratsimp(%)$
normalizeMatrix(%)$
is(% = ident(8));

Result

(%i303) /* degeneracy to Bx=By=Bz=0 eigensystem */
quasiLinearDerivWrtSym_zeroBx_eigenvectors_reduced$
subst([Bx=0,By=0,Bz=0,cf=c],%)$
is(normalizeMatrix(%)=normalizeMatrix(quasiLinearDerivWrtSym_zeroBxyz_eigenvectors_reduced));

Result

(%i306) /* DEGENERACY BY=BZ=0

... therefore cStar^2 = (c^2 + ca^2)/2
... sqrt(cStar^4 - c^2 ca^2) = sqrt(((c^2 + ca^2)/2)^2- c^2 ca^2) = sqrt(c^2/4 + c^2 ca^2 / 2 + ca^2/4 - c^2 ca^2) = sqrt(c^2/4 - c^2 ca^2 / 2 + ca^2/4)
      = sqrt(((c^2 - ca^2)/2)^2) = |c^2 - ca^2|/2
if c^2 > ca^2 then ...
    cf = cStar^2 + sqrt(cStar^4 - c^2 ca^2) = (c^2 - ca^2)/2 + (c^2 - ca^2)/2 = c^2 - ca^2
    cs = cStar^2 - sqrt(cStar^4 - c^2 ca^2) = (c^2 - ca^2)/2 - (c^2 - ca^2)/2 = 0
if c^2 < ca^2 then ... |c^2 - ca^2|/2 = -(c^2 - ca^2)/2
    cf = cStar^2 + sqrt(cStar^4 - c^2 ca^2) = -(c^2 - ca^2)/2 + (c^2 - ca^2)/2 = 0
    cs = cStar^2 - sqrt(cStar^4 - c^2 ca^2) = -(c^2 - ca^2)/2 - (c^2 - ca^2)/2 = ca^2 - c^2
...but I'm pretty sure c^2 > ca^2 can be proved ... (and it is proved in Trangenstein)*/"";

Result

(%i307) quasiLinearDerivWrtSym$
subst([By=0,Bz=0],%)$
quasiLinearDerivWrtSym_zeroByz : %;

Result

(%i310) quasiLinearDerivWrtSym_zeroByz_eigenvector_results : ev(similaritytransform(quasiLinearDerivWrtSym_zeroByz), hermitianmatrix=true);

Result

(%i311) block([A,B],
    [A,B] : makeEigenMatrices(quasiLinearDerivWrtSym_zeroByz_eigenvector_results),
    B : B . diag([sqrt(2),sqrt(2),-sqrt(2),-sqrt(2),sqrt(2),sqrt(2),1,1]),
    [quasiLinearDerivWrtSym_zeroByz_maximaEigenvalues_reduced, quasiLinearDerivWrtSym_zeroByz_maximaEigenvectors_reduced] : [expand(A),B]);

Result

(%i312) permuteEigenvalues_zeroByz : matrix(
[0,0,0,0,1,0,0,0],
[1,0,0,0,0,0,0,0],
[0,1,0,0,0,0,0,0],
[0,0,0,0,0,0,1,0],
[0,0,0,0,0,0,0,1],
[0,0,1,0,0,0,0,0],
[0,0,0,1,0,0,0,0],
[0,0,0,0,0,1,0,0]
);

Result

(%i313) permuteEigenvalues_zeroByz . quasiLinearDerivWrtSym_zeroByz_maximaEigenvalues_reduced . invert(permuteEigenvalues_zeroByz)$
quasiLinearDerivWrtSym_zeroByz_eigenvalues_reduced : %;

Result

(%i315) quasiLinearDerivWrtSym_zeroByz_maximaEigenvectors_reduced . invert(permuteEigenvalues_zeroByz)$
quasiLinearDerivWrtSym_zeroByz_eigenvectors_reduced : %;

Result

(%i317) /* eigenvector orthogonality */
transpose(quasiLinearDerivWrtSym_zeroByz_eigenvectors_reduced) . quasiLinearDerivWrtSym_zeroByz_eigenvectors_reduced$
normalizeMatrix(%)$
is(%=ident(8));

Result

Created with wxMaxima.