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?
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). $$
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.
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.
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}.$$
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)
result = .... Also, to whoever it may help: the expm is from scipy.linalg.
– anderstood
Jan 29 '19 at 14:13
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.
The other answers describe nicely how to compute the integral $$ \int_0^T e^{At}dt$$ but the conditions under which this integral is invertible have not been fully described. The matrix is in fact invertible if an only if, for each non-zero eigenvalues of $A$, let us call it $\lambda$, the condition $$e^{\lambda T} \neq 1$$ holds. This is equivalent to saying that for every eigenvalue $\lambda$ of $A$ we must have $$\lambda \neq 2 \pi k i/T$$ for all non-zero integers $k$.
Explanation:
This does in fact follow from a theorem by Kalman-Ho-Narendra, but not quite in the way obareey describes in another answer. I will repeat some of the steps that show the above from the proof of Theorem 12 in "Controllability of Linear Dynamical Systems" written by R. E. Kalman, Y. C. Ho and K. S. Narendra, published in the journal Contributions to Differential Equations, volume 1 (1963), pp. 189-213.
Let $F$ be the Jordan form of the matrix $A$, so that $A = U^{-1}FU$ for some full-rank matrix $U$. From properties of the matrix exponential, we have $$\int_0^T e^{At} dt = U^{-1} \int_0^T e^{Ft}dt U\\ \implies det\left(\int_0^T e^{At} dt\right) = det(U^{-1})det\left(\int_0^T e^{Ft}dt\right)det(U)$$ Since $U$ is full rank, we have that $\int_0^T e^{At}dt$ is invertible if and only if $\int_0^T e^{Ft}dt$ is also full rank (because then the original integral has a non-zero determinant). Let $F$ consist of $s$ Jordan blocks so that $$F = \begin{bmatrix} J_1(\lambda_1) & \cdots & 0 \\ \vdots & \ddots & \vdots \\ 0 & \cdots & J_s(\lambda_s) \end{bmatrix} \\ \implies e^{Ft} = \begin{bmatrix} e^{J_1(\lambda_1)t} & \cdots & 0 \\ \vdots & \ddots & \vdots \\ 0 & \cdots & e^{J_s(\lambda_s)t} \end{bmatrix},$$ where we used properties of the matrix exponential when applied to block-diagonal matrices. Now, $\int_0^T e^{Ft}dt$ will also be block-diagonal, and full rank if and only if each of its diagonal blocks is full-rank. Note that $$e^{J_i(\lambda_i)t} = e^{\lambda_it} \begin{bmatrix} 1 & t & t^2 & \cdots & \frac{t^{r_i-1}}{(r_i-1)!} \\ 0 & 1 & t & \cdots & \frac{t^{r_i-2}}{(r_i-2)!} \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & 0 & \cdots & 1 \end{bmatrix}$$ where $r_i$ is the size of the Jordan block $J_i(\lambda_i)$. I don't derive this here, but it's easy to find the derivation when searching for the matrix exponential of a Jordan matrix. Now, if we define $$b_k = \int_0^T e^{\lambda_it} \frac{t^k}{k!}dt\\ \implies \int_0^T e^{J_i(\lambda_i)t}dt = \begin{bmatrix} b_0 & b_1 & b_2 & \cdots & b_{r_i-1} \\0 & b_0 & b_1 & \cdots & b_{r_i-2} \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & 0 & \cdots & b_0 \end{bmatrix}$$ where $$b_0 = \int_0^T e^{\lambda_i t}dt = \begin{cases} \frac{e^{\lambda_iT}-1}{\lambda_i}, \lambda_i \neq 0 \\ T, \lambda_i = 0\end{cases}.$$ Now, from properties of triangular matrices, we know that all eigenvalues of $e^{J_i(\lambda_i)}$ are equal to $b_0$. We can apply this process for every diagonal block of $\int_0^T e^{Ft}dt$, and we get that a block is not full rank if and only if $e^{\lambda_iT}=1$. That is, $\int_0^T e^{Ft}dt$ and by extension $\int_0^T e^{At} dt$ are full rank if and only if every eigenvalue $\lambda_i$ of $A$ satisfies $\lambda_i \neq 2\pi k i/T$, proving the above claims.
Sorry for the lengthy explanation, let me know if anything isn't clear!
