Return to Revisions

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$]

• For $A$ and $B$ two matrices $\left(A.B\right)^{-1}=B^{-1}.A^{-1}$

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 (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] *)