1
$\begingroup$

I have a matrix that has a variable number of columns and rows (but it is always square). The matrix is:

$$ \Gamma = \begin{bmatrix} \Theta & 0 & 0 & 0 & T & 0 & R_j \\ u \Theta & \rho & 0 & 0 & u T & 0 & R_j \\ v \Theta & 0 & \rho & 0 & v T & 0 & R_j \\ w \Theta & 0 & 0 & \rho & w T & 0 & R_j \\ H \Theta -1 & u \rho & v \rho & w \rho & \Omega \rho +H T & \frac{5 \rho }{3} & R_j \\ k \Theta & 0 & 0 & 0 & k T & \rho & k R_j \\ Y_i \Theta & 0 & 0 & 0 & Y_i T & 0 & Y_i R_j+\rho \delta_{ij} \\ \end{bmatrix} $$

where the indices $i,j$ are the row and column indices respectively. For example, if $i,j={1,2}$ the matrix would look like:

$$ \Gamma = \begin{bmatrix} \Theta & 0 & 0 & 0 & T & 0 & R_1 & R_2 \\ u \Theta & \rho & 0 & 0 & u T & 0 & R_1 & R_2 \\ v \Theta & 0 & \rho & 0 & v T & 0 & R_1 & R_2 \\ w \Theta & 0 & 0 & \rho & w T & 0 & R_1 & R_2 \\ H \Theta -1 & u \rho & v \rho & w \rho & \Omega \rho +H T & \frac{5 \rho }{3} & R_1 & R_2 \\ k \Theta & 0 & 0 & 0 & k T & \rho & k R_1 & k R_2 \\ Y_1 \Theta & 0 & 0 & 0 & Y_1 T & 0 & Y_1 R_1+\rho & Y_1 R_2 \\ Y_2 \Theta & 0 & 0 & 0 & Y_2 T & 0 & Y_2 R_1 & Y_2 R_2 + \rho \\ \end{bmatrix} $$

The values for $i$ and $j$ can be arbitrary.

When I put this in as is (just treating the subscripts as notation), I get terms divided by the Kronecker delta only which is impossible. Is there a way to invert this and keep it general in terms of $i,j$? The code to enter the matrix is:

G = 
  {
    {Θ, 0, 0, 0, T, 0, Rj},
    {Θ*u, ρ, 0, 0, T*u, 0, Rj},
    {Θ*v, 0, ρ, 0, T*v, 0, Rj},
    {Θ*w, 0, 0, ρ, T*w, 0, Rj},
    {Θ*H - 1, ρ*u, ρ*v, ρ*w, T*H + ρ*Ω, 5/3*ρ, Rj},
    {Θ*k, 0, 0, 0, T*k, ρ, Rj*k},
    {Θ*Yi, 0, 0, 0, T*Yi, 0, Rj*Yi + ρ*δ}
  };
Ginv = FullSimplify[Inverse[G]];
$\endgroup$
  • $\begingroup$ Can you explain what is wrong with this? Ginv.G gives the right thing. Also, what does "terms divided by the Kronecker delta only" mean, and why is it impossible? $\endgroup$ – acl Jul 10 '14 at 22:51
  • $\begingroup$ Mathematica does not know what these symbols mean mathematically. there are just symbols. If you want actual dirac delta, you should use DiracDelta[] btw, you might want to add space here "[Rho][Delta]" to make it "[Rho] [Delta]" just in case. $\endgroup$ – Nasser Jul 10 '14 at 22:53
  • $\begingroup$ @acl If I look at the Ginv, there are terms like $Y_i/\delta_{ij}$ which will be $Y_i/0$ for all but 1 column, which is ill-defined. $\endgroup$ – tpg2114 Jul 10 '14 at 22:54
  • $\begingroup$ @Nasser But how can I tell Mathematica that $i,j$ are for columns and rows? Without that information, the DiracDelta[] wouldn't know what to do either right? $\endgroup$ – tpg2114 Jul 10 '14 at 22:55
  • 1
    $\begingroup$ I do not understand what you have there, so can't help. I do not know what Y_i supposed to mean. And do not know what R_j suppose to mean. Where did these come from? What is "i" there? what is "j" ? your question is fully described well to answer it (for me at least) $\endgroup$ – Nasser Jul 10 '14 at 23:03
2
$\begingroup$

This seems to work although I have not tested it thoroughly. Also, it's a rather brute force approach, so I hope someone else will post a more elegant solution.

g[0] = 
  {{Θ, 0, 0, 0, T, 0}, 
   {u Θ, ρ, 0, 0, T u, 0}, 
   {v Θ, 0, ρ, 0, T v, 0}, 
   {w Θ, 0, 0, ρ, T w, 0}, 
   {-1 + H Θ, u ρ, v ρ, w ρ, H T + ρ Ω, (5 ρ)/3}, 
   {k Θ, 0,0, 0, k T, ρ}}

upperRight[n_] := Transpose@Table[Join[ConstantArray[R[j], 5], {k R[j]}], {j, n}]

lower[n_] := 
  Table[Join[{Θ*Y[i], 0, 0, 0, T*Y[i], 0}, Table[Y[i] R[j], {j, n}]], {i, n}] /. 
    p : R[x_] Y[x_] :> p + ρ

g[n_] := Join[MapThread[Join, {g[0], upperRight[n]}, 1], lower[n]]

g2 = g[2]
{{Θ, 0, 0, 0, T, 0, R[1], R[2]}, 
  {u Θ, ρ, 0, 0, T u, 0, R[1], R[2]}, 
  {v Θ, 0, ρ, 0, T v, 0, R[1], R[2]}, 
  {w Θ, 0, 0, ρ, T w, 0, R[1], R[2]}, 
  {-1 + H Θ, u ρ, v ρ, w ρ, H T + ρ Ω, (5 ρ)/3, R[1], R[2]}, 
  {k Θ, 0, 0, 0, k T, ρ, k R[1], k R[2]}, 
  {Θ Y[1], 0, 0, 0, T Y[1], 0, ρ + R[1] Y[1], R[2] Y[1]}, 
  {Θ Y[2], 0, 0, 0, T Y[2], 0, R[1] Y[2], ρ + R[2] Y[2]}}

The above is invertible although the result is messy.

Short[Inverse[g2], 4]

matrix

$\endgroup$
1
$\begingroup$

As it turns out, by employing block matrix inversion, one can derive a general formula for the inverse of $\mathbf \Gamma$. I'll leave the rigorous algebra/justification for other people, but here is how to apply block inversion:

emat = {{Θ, 0, 0, 0, T, 0},
        {Θ u, ρ, 0, 0, T u, 0},
        {Θ v, 0, ρ, 0, T v, 0},
        {Θ w, 0, 0, ρ, T w, 0},
        {Θ H - 1, ρ u, ρ v, ρ w, T H + ρ Ω, 5 ρ/3},
        {Θ k, 0, 0, 0, T k, ρ}};

n = 7; (* size of lower right block, change as needed *)

(* upper right *)
ff = DiagonalMatrix[{1, 1, 1, 1, 1, k}].
     Transpose[Table[ConstantArray[Subscript[R, j], 6], {j, n}]];

(* lower left *)
gg = Table[{Subscript[Y, i] Θ, 0, 0, 0, Subscript[Y, i] T, 0}, {i, n}];

(* lower right *)
hh = Outer[Times, Table[Subscript[Y, i], {i, n}], Table[Subscript[R, j], {j, n}]] +
     ρ IdentityMatrix[n];

Generate the full matrix $\mathbf \Gamma$:

gmat = ArrayFlatten[{{emat, ff}, {gg, hh}}];

Now, assemble the inverse in blocks:

einv = Simplify[Inverse[emat]]; (* invert upper left block *)

vec = Simplify[Total[einv.DiagonalMatrix[{1, 1, 1, 1, 1, k}], {2}]];

dp = Table[Subscript[R, j], {j, n}].Table[Subscript[Y, j], {j, n}];

Build $\mathbf \Gamma^{-1}$:

ginv = ArrayFlatten[{{einv + Transpose[PadRight[{vec dp/ρ}, {6, 6}]],
                      -Outer[Times, vec, Table[Subscript[R, j], {j, n}]]/ρ},
                     {-Table[PadRight[{Subscript[Y, i]}, 6], {i, n}]/ρ,
                      IdentityMatrix[n]/ρ}}];

Check:

gmat.ginv == IdentityMatrix[6 + n] // Simplify
   True

In symbolic format, the inverse can be expressed as

$$\mathbf\Gamma^{-1}=\small \begin{pmatrix} \frac{3 H T-5 k T-3 \left(u^2+v^2+w^2\right)T+3\rho\Omega}{3(T+\Theta\rho\Omega)}+\frac{\left(T\left(H-u^2+u-v^2+v-w^2+w-1\right)+\rho\Omega\right)\sum_k R_k Y_k}{\rho(\Theta\rho\Omega+T)}&\frac{T u}{T+\Theta\rho\Omega} &\frac{T v}{T+\Theta\rho\Omega}&\frac{T w}{T+\Theta\rho \Omega}&-\frac{T}{T+\Theta\rho\Omega}&\frac{5T}{3(T+\Theta\rho\Omega)}&-\frac{T\left(H-u^2+u-v^2+v-w^2+w-1\right)+\rho\Omega}{\rho(\Theta\rho\Omega+T)}R_j\\ \frac{(1-u)\sum_k R_k Y_k}{\rho^2}-\frac{u}{\rho}&\frac{1}{\rho} & 0 & 0 & 0 & 0 &\frac{u-1}{\rho^2}R_j\\ \frac{(1-v)\sum_k R_k Y_k}{\rho^2}-\frac{v}{\rho} & 0 & \frac{1}{\rho} & 0 & 0 & 0 &\frac{v-1}{\rho^2}R_j\\ \frac{(1-w)\sum_k R_k Y_k}{\rho^2}-\frac{w}{\rho} & 0 & 0 & \frac{1}{\rho} & 0 & 0 &\frac{w-1}{\rho^2}R_j\\ \frac{-3 H\Theta+5 k\Theta+3\left(u^2+v^2+w^2\right)\Theta+3}{3(T+\Theta\rho\Omega)}+\frac{\left(1-\Theta\left(H-u^2+u-v^2+v-w^2+w-1\right)\right)\sum_k R_k Y_k}{\rho(\Theta\rho\Omega +T)}& -\frac{u\Theta}{T+\Theta\rho\Omega}& -\frac{v\Theta}{T+\Theta\rho\Omega} & -\frac{w\Theta}{T+\Theta\rho\Omega} & \frac{\Theta}{T+\Theta\rho\Omega}& -\frac{5\Theta}{3 (T+\Theta\rho\Omega)}&-\frac{1-\Theta\left(H-u^2+u-v^2+v-w^2+w-1\right)}{\rho(\Theta\rho\Omega +T)}R_j\\ -\frac{k}{\rho} & 0 & 0 & 0 & 0 & \frac{1}{\rho} & 0 \\ -\frac{Y_i}{\rho} & 0 & 0 & 0 & 0 & 0 & \frac{\delta_{ij}}{\rho} \end{pmatrix} $$


The block inversion worked as well as it did because of fortuitous simplifications. For instance, here are some useful identities I determined:

einv.ff == Outer[Times, vec, Table[Subscript[R, j], {j, n}]] // Simplify
   True

gg.einv == PadRight[Table[{Subscript[Y, i]}, {i, n}], {n, 6}] // Simplify
   True

Inverse[hh - gg.einv.ff] == IdentityMatrix[n]/ρ // Simplify
   True

It was by using these identities that I was able to construct the explicit expression for $\mathbf \Gamma^{-1}$.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.