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?

Rodrigo de Azevedo
  • 18,977
  • 5
  • 36
  • 95
  • 431
  • 1
  • 4
  • 4
  • 1
    you can represent matrix exponential by power series – dato datuashvili Jan 31 '14 at 09:52
  • 1
    http://en.wikipedia.org/wiki/Matrix_exponential – dato datuashvili Jan 31 '14 at 09:52
  • 1
    Spectral theorem says that if you take analytical function $f$, apply it to the matrix $A$ (its eigenvaules has to be in the domain of analaticity of $f$), then eigenvalues of $f(A)$ are $f(\text{eigenvalues of } A)$. Exponent is never zero, hence $\exp(A)$ is always invertible, moreover, $\exp(A)^{-1}=\exp(-A)$. – TZakrevskiy Jan 31 '14 at 10:05

5 Answers5


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). $$

Yiorgos S. Smyrlis
  • 78,494
  • 15
  • 113
  • 210

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.

  • 5,404
  • 1
  • 15
  • 32

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 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)

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.

  • 3,414
  • 13
  • 32