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?

  • 1
    $\begingroup$ you can represent matrix exponential by power series $\endgroup$ Commented Jan 31, 2014 at 9:52
  • 1
    $\begingroup$ en.wikipedia.org/wiki/Matrix_exponential $\endgroup$ Commented Jan 31, 2014 at 9:52
  • 1
    $\begingroup$ 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)$. $\endgroup$ Commented Jan 31, 2014 at 10:05

5 Answers 5


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

  • $\begingroup$ What we will do if we use $A(t)$ instead of $tA$ $\endgroup$
    – Nirvana
    Commented Oct 6, 2014 at 8:06
  • $\begingroup$ There is nothing wrong with that. $\endgroup$ Commented Oct 6, 2014 at 10:22
  • $\begingroup$ @ Yiorgos S. Smyrlis So you are saying even if $A(t)$ and $A^{'}(t)$ are not commutative,this will work $\endgroup$
    – Nirvana
    Commented Oct 6, 2014 at 10:28
  • $\begingroup$ Here A is not a constant matrix. It is a variable matrix $\endgroup$
    – Nirvana
    Commented Oct 6, 2014 at 10:28
  • 1
    $\begingroup$ If $A$ is a psd matrix, can we say $\int_0^T e^{tA}\,dt = (e^{TA}-I)A^\dagger$, where $A^\dagger$ is the Moore-Penrose pseudo-inverse of $A$ ? $\endgroup$
    – dohmatob
    Commented Jun 1, 2022 at 18:15

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.

  • $\begingroup$ Thank you. This answer looks very promising. Btw, I am looking at this matrix exactly because of what you said about the discretization! $\endgroup$
    – samsa44
    Commented Feb 24, 2014 at 11:19
  • $\begingroup$ I just want to know what's the reference book or paper for Yiorgios or Obareey? $\endgroup$ Commented Dec 7, 2020 at 19:57
  • $\begingroup$ You can finish the answer by writing the integral equals $A^{-1}(e^{AT}-I)$ $\endgroup$ Commented Sep 10, 2023 at 12:43
  • $\begingroup$ @ТymaGaidash٠ yes, this works when $A$ is invertible. The answer does not make this assumption. $\endgroup$
    – obareey
    Commented Sep 11, 2023 at 7:17

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

  • $\begingroup$ Very nice implementation and explanation. $\endgroup$
    – Clay
    Commented Jul 7, 2021 at 19:05
  • $\begingroup$ +1 Very nice explanation. $\endgroup$
    – greg
    Commented Sep 22, 2022 at 22:21
  • $\begingroup$ What would happen if you multiplied out the vectors? $\endgroup$ Commented Sep 10, 2023 at 12:40
  • $\begingroup$ Isn't the solution to the non-homogeneous ODE: $X(t) = \mathrm{e}^{At}X(0) + \int_0^t\mathrm{e}^{A(t-\tau)}\mathrm{d}\tau B$? (see difference in exponent of integrand...) $\endgroup$
    – jmd
    Commented Jun 8 at 22:59

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)
  • $\begingroup$ I get a shape error in result = ... . Also, to whoever it may help: the expm is from scipy.linalg. $\endgroup$
    – anderstood
    Commented Jan 29, 2019 at 14:13

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.


You must log in to answer this question.

Not the answer you're looking for? Browse other questions tagged .