If we have a continuous time state-equation,
$$ \dot{x}(t) = A x(t) + B u(t)$$
where $A \in \mathbb{R}^{n \times n}, x \in \mathbb{R}^{n\times1}, B \in \mathbb{R}^{n \times m}, u \in \mathbb{R}^{m \times 1}$, the analytical solution is obtained as
$$x(t) = e^{A t} x(0) + \int_0^t e^{A(t-\tau)}B u(\tau)\, d\tau$$
Converting this to discrete-time, with sampling period $T_s$, and assuming that all eigenvalues of $A$ are non-zero, real & distinct and $u$ remains constant within each sample interval, $k$ to $k+1$, we get the following difference equation,
$$x[k+1] = e^{A T_s} x[k] + \left[\int_0^{T_s} e^{A\tau}B \, d\tau\right] u[k]$$
which can be written as
$$x[k+1] = A_d x[k] + B_d u[k]$$
When the original matrix $A$ is invertible, the integral term corresponding to $B_d$ can be written as
$$B_d = A^{-1}(A_d - I)$$
The question is, What about the case when A is non-invertible? In particular, if
$$A=\begin{bmatrix}a_1 & 0 & 0 \\ 0 & a_2 & 0 \\ 0 & 0 & 0 \end{bmatrix}$$ and $$B=\begin{bmatrix}b_1 \\ b_2 \\ b_3 \end{bmatrix}$$
then one of the eigenvalues of $A$ is zero. In this case, how do I obtain the corresponding $B_d$?