This answer is an attempt to provide a solution for a hypothetical user that does not know the general formula beforehand and is unaware of any packages.
The answer below considers that the user is aware of the following:
If $f$ is an expandable function that can be considered on the real line (for example $f(x)=(1+x)^{-1}$ or $f(x)=\exp(x)$), then, for a matrix variable $A$ and a real number valued variable $a$, $f(A)$=
Replace
[$f(a),a \rightarrow A$]. This point allows us to use Series for real valued variables.For $A$ and $B$ two matrices $\left(A.B\right)^{-1}=B^{-1}.A^{-1}$. This point is used to place the initial expression into the form $f(x)$ where $f(x)=(1+x)^{-1}$ and $ x=A^{-1}.B $
Then if we define an identity tensor for the dot product:
Dot[s___, id, g___] ^:= Dot[s, g]
The formula may be obtained with the following code that attempts to do simple replacements sequentially (var
below is a formal wrapper that formally replaces a matrix with a real valued variable) :
Fold[ReplaceAll,
Inverse[A + B] ,
{
Inverse[w_] :> Inverse[A . TensorExpand[Inverse@A . w]] ,
Inverse[a_ . b_] :> Inverse[b] . Inverse[a],
MatrixPower[_, 0] -> id,
Inverse[id + s_] :> ReplaceAll[
Normal@Series[(1 + epsilon*var[s])^(-1),
{epsilon, 0, 3}] /. epsilon -> 1
, {
var[a_]^m_:>Dot@@ConstantArray[a,m]
, var->Identity
, 1 -> id
}
]
}
]//TensorExpand
out: (* -MatrixPower[A, -1] . B . MatrixPower[A, -1] + MatrixPower[A, -1] . B . MatrixPower[A, -1] . B . MatrixPower[A, -1] - MatrixPower[A, -1] . B . MatrixPower[A, -1] . B . MatrixPower[A, -1] . B . MatrixPower[A, -1] + MatrixPower[A, -1] *)