(%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 */"";

 (%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);

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

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

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

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

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

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

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

 (%i16) vSq : vs.vs;

 (%i17) BSq : Bs.Bs;

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

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

 (%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);

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

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

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

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

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

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

 (%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);

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

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

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

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

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

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

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

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

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

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

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

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

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

 (%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);

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

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

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

 (%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 : %;

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

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

 (%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 : %;

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

 (%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 : %;

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

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

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

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

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

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

 (%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 : %;

 (%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);

 (%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 : %;

 (%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 : %;

 (%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 : %;

 (%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));

 (%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;

 (%i134) dPrim_dSym : invert(dSym_dPrim);

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

 (%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 : %;

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

 (%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 */"";

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

 (%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 : %;

 (%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 : %;

 (%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));

 (%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));

 (%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 */"";

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

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

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

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

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

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

 (%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 */"";

 (%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] );

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

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

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

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

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

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

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

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

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

 (%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 :%;

 (%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],%) */ %;

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

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

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

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

 (%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]);

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

 (%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]);

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

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

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

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

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

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

 (%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 : %;

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

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

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

 (%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] );

 (%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]);

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

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

 (%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 : %;

 (%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.*/"";

 (%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));

 (%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));

 (%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)*/"";

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

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

 (%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]);

 (%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] );

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

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

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

Created with wxMaxima.