Nonhomogenous matrix differential equation

Question

In this paper, I've stumbled unpon the solution they give for an non-homogeneous matrix differential equation.

Basically, the equation is:

$\frac{\partial\boldsymbol{\psi(t)}}{\partial t} = \boldsymbol{A\psi}(t) + \boldsymbol{B}$

Here, A is a constant matrix and B a constant vector, between $t$ and $t+\Delta t$.

Now, they give in equation (12) that the general solution for this equation is:

$\boldsymbol{\psi} (t+\Delta t) = e^{\boldsymbol{A}\Delta t} \boldsymbol{\psi} (t) + \left(e^{\boldsymbol{A}\Delta t} - \boldsymbol{I} \right) \boldsymbol{A^{-1}} \boldsymbol{B}$

I have been unable to find the derivation for this in the literature. Where does the inverse matrix A come from?

Do you know of a document explaining the derivation of this type of matrix differential equation with a constant second term (I was able to find the easier case where B = 0, and the more complex case where B = U(t), but no "in-between").

My best guess is that we have somehow:

$\boldsymbol{\psi} (t+\Delta t) = e^{\boldsymbol{A}\Delta t} \boldsymbol{\psi} (t) + e^{\boldsymbol{A}\Delta t} *\boldsymbol{B} = e^{\boldsymbol{A}\Delta t} \boldsymbol{\psi} (t) + \left(e^{\boldsymbol{A}\Delta t} - \boldsymbol{I} \right) \boldsymbol{A^{-1}} \boldsymbol{B}$

And that:

$e^{\boldsymbol{A}\Delta t} *\boldsymbol{B}$ is the convolution product of the exponential of the matrix with (in this case) a constant vector, and that this convolution product can be written $\left(e^{\boldsymbol{A}\Delta t} - \boldsymbol{I} \right) \boldsymbol{A^{-1}} \boldsymbol{B}$, though I'm not sure how and haven't been able to find information on this (maybe because I'm not using exactly the right terminology)

Answer 1

The derivation of the solution to

$\dot \psi(t) = A\psi(t) + B \tag 1$

given that $\psi(t)$ takes the value $\psi(t_0)$ at $t = t_0$, is in fact very well-known. I'll give a reference later, but first I'll try to present a quick run-down of how the formula for $\psi(t)$ may be established. (1) may be written in the form

$\dot \psi(t) - A\psi(t) = B, \tag 2$

and multiply through by the matrix exponential $e^{-A(t - t_0)}$:

$e^{-A(t - t_0)} \dot \psi - A e^{-A(t - t_0)} \psi = e^{-A(t - t_0)}B, \tag 3$

and observe that the left-had side now has an anti-deriviative; that is,

$e^{-A(t - t_0)} \dot \psi - A e^{-A(t - t_0)} \psi = \dfrac{d}{dt}(e^{-A(t - t_0)}\psi) = e^{-A(t - t_0)}B; \tag 4$

this equation may be integrated 'twixt $t_0$ and $t$:

$e^{-A(t - t_0)}\psi(t) - \psi(t_0) = e^{-A(t - t_0)}\psi(t) - e^{-A(t_0 - t_0)}\psi(0) = \displaystyle \int_{t_0}^t \dfrac{d}{ds}(e^{-A(s - t_0)}\psi(s)) \ ds = \int_{t_0}^t e^{-As}B \; ds = \left (\int_{t_0}^t e^{-As} \; ds \right ) B; \tag 5$

it is essy to verify via direct differentiation that

$\displaystyle \int_{t_0}^t e^{-As} \; ds = A^{-1}(I - e^{=A(t- t_0)}); \tag 6$

I leave the details to the reader. Substituting (6) into (5) yields

$e^{-A(t - t_0)}\psi(t) - \psi(t_0) = A^{-1}(I - e^{=A(t - t_0)})B, \tag 7$

or

$e^{-A(t - t_0)}\psi(t) = \psi(t_0) + A^{-1}(I - e^{=A(t - t_0)})B, \tag 8$

which may now be multiplied through by $e^{A(t - t_0)}$:

$\psi(t) = e^{A(t - t_0)}\psi(t_0) + A^{-1}(e^{A(t - t_0)} - I)B; \tag 9$

Check:

from (9),

$\dot \psi(t) = Ae^{A(t - t_0)}\psi(t_0) + A^{-1}Ae^{A(t - t_0)}B = Ae^{A(t - t_0)}\psi(t_0) + AA^{-1}e^{A(t - t_0)}B = A(e^{A(t - t_0)}\psi(t_0) + A^{-1}e^{A(t - t_0)}B)$ $= A(e^{A(t - t_0)}\psi(t_0) + A^{-1}(e^{A(t - t_0)} - I)B) + B = A\psi(t) + B. \tag{10}$

For more information, see this wikipedia page.

Answer 2

Note first that we do have to assume that $\boldsymbol{A}$ is nonsingular, because that's the only way the formula can be well-defined as written. We have the matrix differential equation $$ \frac{\partial\boldsymbol{\psi(t)}}{\partial t} = \boldsymbol{A\psi}(t) + \boldsymbol{B}$$ We can use homogeneous coordinates, which are a somewhat clever way to make this affine transformation of $\boldsymbol{\psi}$ linear in some $\boldsymbol{\Psi}$. Writing $$ \boldsymbol{\Psi}(t) = \begin{bmatrix} \boldsymbol{\psi(t)} \\ \boldsymbol{I}\end{bmatrix} $$ we have $$\frac{\partial \boldsymbol{\Psi(t)}}{\partial t} = \frac{\partial}{\partial t}\begin{bmatrix}\boldsymbol{\psi(t)} \\ \boldsymbol{I} \end{bmatrix} = \begin{bmatrix} \frac{\partial \boldsymbol{\psi(t)}}{\partial t} \\ \boldsymbol{0} \end{bmatrix} = \begin{bmatrix} \boldsymbol{A\psi}(t) + \boldsymbol{B}\\ \boldsymbol{0} \end{bmatrix} = \underbrace{\begin{bmatrix} \boldsymbol{A} & \boldsymbol B \\ \boldsymbol{0} & \boldsymbol{0} \end{bmatrix}}_{\boldsymbol{\Phi}}\begin{bmatrix}\boldsymbol{\psi(t)} \\ \boldsymbol{I} \end{bmatrix} = \boldsymbol{\Phi \Psi(t)}$$ which is the sort of homogeneous equation we know how to solve. In particular, it has solution $$ \boldsymbol{\Psi(t + \Delta t)} = \exp\left(\boldsymbol{\Phi}\Delta t \right)\boldsymbol{\Psi}(t) $$ You've said you can solve the case $\boldsymbol{B=0}$, so I'm assuming this solution for homogeneous equations is familiar to you. Now, one can consider the pattern and eventually show by induction that for $k \geq 1$ $$\begin{bmatrix} \boldsymbol{A} & \boldsymbol B \\ \boldsymbol{0} & \boldsymbol{0} \end{bmatrix}^k = \begin{bmatrix} \boldsymbol{A}^k & \boldsymbol A^{k-1}\boldsymbol{B} \\ \boldsymbol{0} & \boldsymbol{0} \end{bmatrix}$$ Using the Taylor series definition of the exponential then, on gets $$\exp\left( \boldsymbol{\Phi}\right) = \boldsymbol{I} + \sum_{k=1}^\infty \left(\frac{1}{k!}\begin{bmatrix} \boldsymbol{A}^k & \boldsymbol A^{k-1}\boldsymbol{B} \\ \boldsymbol{0} & \boldsymbol{0} \end{bmatrix} (\Delta t)^k \right)$$ and most of the blocks of our exponential are easy. The top left block is just the series definition of $\exp(\boldsymbol{A}\Delta t)$. And the bottom two blocks in the sum above are zero, so that in $\exp(\boldsymbol{\Phi} \Delta t)$, the bottom left block is just $\boldsymbol{0}$ and the bottom right block is just $\boldsymbol{I}$. So know we only need to consider the top right block. Using the distributive property for matrix-vector multiplication, we have that the top-right block is just $\boldsymbol{B}$ left multiplied by the infinite matrix sum (let's call it $\boldsymbol{S}$) $$\boldsymbol{S} = \sum_{k=1}^\infty \left(\frac{(\Delta t)^k}{k!}\boldsymbol{A}^{k-1}\right) $$ We have $$ \boldsymbol{SA} = \sum_{k=1}^\infty \left(\frac{(\Delta t)^k}{k!}\boldsymbol{A}^{k}\right) $$ and therefore $$ \boldsymbol{SA} + \boldsymbol{I} = \sum_{k=0}^\infty \left(\frac{(\Delta t)^k}{k!}\boldsymbol{A}^{k}\right) = \exp(\boldsymbol{A}\Delta t)$$ so that, using our assumption $\boldsymbol{A}$ is invertible $$ \boldsymbol{S} = \left( e^{\displaystyle{\boldsymbol{A}\Delta t}} - \boldsymbol{I}\right) \boldsymbol{A}^{-1}$$ and in particular the top-right block of $\exp \left(\boldsymbol{\Phi}\Delta t\right)$ is just $$ \boldsymbol{SB} = \left( e^{\displaystyle{\boldsymbol{A}\Delta t}} - \boldsymbol{I}\right) \boldsymbol{A}^{-1} \boldsymbol{B}$$ This means we have $$ \exp\left( \boldsymbol{\Phi} \Delta t\right) = \begin{bmatrix} \displaystyle{e^{\displaystyle{\boldsymbol{A}\Delta t}}} & \displaystyle{\left( e^{\displaystyle{\boldsymbol{A}\Delta t}} - \boldsymbol{I}\right) \boldsymbol{A}^{-1} \boldsymbol{B}} \\ \boldsymbol{0} & \boldsymbol{I}\end{bmatrix} $$ so that we have from the solution to homogeneous case again that $$ \boldsymbol{\Psi}(t + \Delta t) = \exp (\boldsymbol{\Phi}\Delta t) \boldsymbol{\Psi}(t) = \begin{bmatrix} \displaystyle{e^{\displaystyle{\boldsymbol{A}\Delta t}}} & \displaystyle{\left( e^{\displaystyle{\boldsymbol{A}\Delta t}} - \boldsymbol{I}\right) \boldsymbol{A}^{-1} \boldsymbol{B}} \\ \boldsymbol{0} & \boldsymbol{I}\end{bmatrix} \begin{bmatrix} \boldsymbol{\psi}(t) \\ \boldsymbol{I}\end{bmatrix}$$ so that we have $$\begin{bmatrix} \boldsymbol{\psi}(t+\Delta t) \\ \boldsymbol{I} \end{bmatrix} = \begin{bmatrix} \displaystyle{e^{\displaystyle{\boldsymbol{A}\Delta t}}} \boldsymbol{\psi}(t)+ \displaystyle{\left( e^{\displaystyle{\boldsymbol{A}\Delta t}} - \boldsymbol{I}\right) \boldsymbol{A}^{-1} \boldsymbol{B}}\\ \boldsymbol{I} \end{bmatrix} $$ which, comparing top entries, gives your desired identity $$ \boxed{\boldsymbol{\psi}(t+\Delta t) = \displaystyle{e^{\displaystyle{\boldsymbol{A}\Delta t}}} \boldsymbol{\psi}(t)+ \displaystyle{\left( e^{\displaystyle{\boldsymbol{A}\Delta t}} - \boldsymbol{I}\right) \boldsymbol{A}^{-1} \boldsymbol{B}}}$$

P.S. A lot of the $\boldsymbol{I}$s above are really just $1$s in the case where $\boldsymbol{\psi}$ is a vector. But the use of the identity generalizes to certain cases (in particular, to the one where $\boldsymbol{\psi}$ and $\boldsymbol{B}$ are square matrices, I think, among others), so we write $\boldsymbol{I}$ slightly more generally.