Integral of matrix exponential Let $A$ be an $n \times n$ matrix. Then the solution of the initial value problem
\begin{align*}
  \dot{x}(t) = A x(t), \quad x(0) = x_0
\end{align*}
is given by $x(t) = \mathrm{e}^{At} x_0$.
I am interested in the following matrix
\begin{align*}
   \int_{0}^T \mathrm{e}^{At}\, dt
\end{align*}
for some $T>0$. Can one write down a general solution to this without distinguishing cases (e.g. $A$ nonsingular)?
Is this matrix always invertible?
 A: The way I like to do it is based on the following observation: let
$$ \bar{A} := 
     \begin{bmatrix} 
        A & B \\
        0 & 0
     \end{bmatrix}, $$
where $0$ is the zero matrix (dimensions s.t. $\bar{A}$ is square). Then,
$$ \mathrm{e}^{\bar{A}t} = \begin{bmatrix} \mathrm{e}^{At} & \int_0^t\mathrm{e}^{A\tau}\mathrm{d}\tau B \\ 0 & I \end{bmatrix}. $$
Hence, for the integral you can just build this block matrix with $B=I$, compute the matrix exponential of it, then extract the top right block. In a more "closed" form:
$$ \int_0^t\mathrm{e}^{A\tau}\mathrm{d}\tau B = \begin{bmatrix}I & 0\end{bmatrix}\mathrm{e}^{\bar{A}t}\begin{bmatrix}0 \\ I\end{bmatrix}. $$
The advantage of this method with respect to using matrix inversion and/or Jordan form is that this method is numerically stable even when $A$ is singular (or close to singular). The disadvantage, of course, is that it takes a 4x bigger matrix as an input.
Why it works
It follows from this observation: if you have the non-homogeneous ODE
$$ \dot{X}(t) = AX(t) + B, $$
its solution is
$$ X(t) = \mathrm{e}^{At}X(0) +  \int_0^t\mathrm{e}^{A\tau}\mathrm{d}\tau B. $$
Define the auxiliary variable $U(t)$, which is constant (i.e., $U(t) = U(0)$ for all positive $t$). Then $\dot{U}(t) = 0$ and we have the system of ODEs
\begin{align*}
\dot{X}(t) &= AX(t) + BU(t), \\
\dot{U}(t) &= 0,
\end{align*}
which is a homogeneous ODE on the augmented variable $\begin{bmatrix} X(t) \\ U(t) \end{bmatrix}$. Therefore, we have
$$\begin{bmatrix} \dot{X}(t) \\ \dot{U}(t) \end{bmatrix} = \begin{bmatrix} 
        A & B \\
        0 & 0
     \end{bmatrix}\begin{bmatrix} {X}(t) \\ {U}(t) \end{bmatrix} = \bar{A}\begin{bmatrix} {X}(t) \\ {U}(t) \end{bmatrix}$$
whose solution is
$$\begin{bmatrix} {X}(t) \\ {U}(t) \end{bmatrix} = \mathrm{e}^{\bar{A}t}\begin{bmatrix} {X}(0) \\ {U}(0) \end{bmatrix},$$
but also,
$$\begin{bmatrix} {X}(t) \\ {U}(t) \end{bmatrix} = \begin{bmatrix} \mathrm{e}^{At}X(0) +  \int_0^t\mathrm{e}^{A\tau}\mathrm{d}\tau BU(0) \\ U(0) \end{bmatrix} = \begin{bmatrix} \mathrm{e}^{At} &  \int_0^t\mathrm{e}^{A\tau}\mathrm{d}\tau B \\ 0 & I \end{bmatrix}\begin{bmatrix} X(0) \\ U(0) \end{bmatrix}.$$
A: A Python numerical answer
It is surprisingly difficult to find a proper python package for numerical integration of matrix. I know it is not what the question want but I cannot find anywhere else to publish this.
Here, I provide a numerical solution to it. Just call the function 
intergral_result = compute_exp_matrix_intergration(A,T)

will be enough 
import numpy as np
def compute_exp_matrix_intergration(A,T,nbins=100):
    f = lambda x: expm(A*x)
    xv = np.linspace(0,T,nbins)
    result = np.apply_along_axis(f,0,xv.reshape(1,-1))
    return np.trapz(result,xv)

A: Case I. If $A$ is nonsingular, then
$$
\int_0^T\mathrm{e}^{tA}\,dt=\big(\mathrm{e}^{TA}-I\big)A^{-1},
$$
where $I$ is the identity matrix.
Case II. If $A$ is singular, then using the Jordan form we can write $A$ as
$$
A=U^{-1}\left(\begin{matrix}B&0\\0&C\end{matrix}\right)U, 
$$
where $C$ is nonsingular, and $B$ is strictly upper triangular. Then
$$
\mathrm{e}^{tA}=U^{-1}\left(\begin{matrix}\mathrm{e}^{tB}&0\\0&\mathrm{e}^{tC}
\end{matrix}\right)U,
$$
and
$$
\int_0^T\mathrm{e}^{tA}\,dt=U^{-1}\left(\begin{matrix}\int_0^T\mathrm{e}^{tB}dt&0\\0&C^{-1}\big(\mathrm{e}^{TC}-I\big)
\end{matrix}\right)U
$$
But $\int_0^T\mathrm{e}^{tB}dt$ may have different expressions. For example if
$$
B_1=\left(\begin{matrix}0&0\\0&0\end{matrix}\right), \quad
B_2=\left(\begin{matrix}0&1\\0&0\end{matrix}\right),
$$
then
$$
\int_0^T\mathrm{e}^{tB_1}dt=\left(\begin{matrix}T&0\\0&T\end{matrix}\right), \quad
\int_0^T\mathrm{e}^{tB_2}dt=\left(\begin{matrix}T&T^2/2\\0&T\end{matrix}\right).
$$
A: The general formula is the power series
$$ \int_0^T e^{At} dt = T \left( I + \frac{AT}{2!} + \frac{(AT)^2}{3!} + \dots + \frac{(AT)^{n-1}}{n!} + \dots \right) $$
Note that also
$$ \left(\int_0^T e^{At} dt \right) A + I = e^{AT} $$
is always satisfied.
A sufficient condition for this matrix to be non-singular is the so-called Kalman-Ho-Narendra Theorem, which states that the matrix $\int_0^T e^{At} dt$ is invertible if
$$ T(\mu - \lambda) \neq 2k \pi i $$
for any nonzero integer $k$, where $\lambda$ and $\mu$ are any pair of eigenvalues of $A$.
Note to the interested: This matrix also comes from the discretization of a continuous linear time invariant system. It can also be said that controllability is preserved under discretization if and only if this matrix has an inverse.
A: Another Python implementation
Here is another implementation in Python, if it can be helpful to anyone... (and since ArtificiallyIntelligence's answer returns an error in my setup).
The value of the integral is integral and the last line verifies the equality $\int_0^T e^{At}dt = A^{-1}(e^{tA}-I)$, provided $A$ is nonsingular (which is generically the case for randomly generated matrices).
import numpy as np
N = 5
t = 1
A = np.random.rand(N,N)
taylor = t*np.array([np.linalg.matrix_power(A*t,k)/np.math.factorial(k+1) for k in range(50)])
integral = taylor.sum(axis = 0)

print np.linalg.norm(integral - np.dot(np.linalg.inv(A),scipy.linalg.expm(t*A)-np.identity(N)))

Note that you should adjust the 50 in taylor = ... to check convergence.
