# Compute state transition matrix for discrete system

Given a discrete-time linear time-invariant (LTI) system like

$$x(k+1) = A x(k)$$

the solution at a given instant $$k$$ can be computed using the state transition matrix $$A^k$$ by

$$x(k) = A^k x(0)$$

For continuous-time LTI systems, the state transition matrix can be computed using the matrix exponential, which gives the continuous state transition matrix with high accuracy.

However, for discrete-time systems, the state transition matrix is not computed with the matrix exponential but with $$A^k$$. In Matlab, for example, I could compute $$A^k$$ directly using a syntax like A^k, which I guess just performs $$k$$ subsequent matrix multiplications. However, I am not sure if this is the best approach if $$k$$ is large.

Question: What is the most efficient way (in terms of accuracy) to compute $$A^k$$ for large $$k$$ numerically (preferable in Matlab)?

• I’m a bit rusty on this, but could you diagonalize and raise that diagonal matrix to the $k$th power? – rb612 Mar 24 at 21:25
• The most efficient depends on the matrix $A$ and what you mean by efficient. Note that $A^{2n} = (A^n)^2$, so some strength reduction can take place as well. If the matrix is unitarily diagonalisable then a numerically nice way is straightforward to implement. – copper.hat Mar 24 at 23:02
• @rb612 can you diagonalize $\begin{bmatrix} 0 & 1 \\ 0 & 0 \end{bmatrix}$? – Mortified Through Math Mar 29 at 16:20

The matrix $$A^k$$ is not what most people refer to as the state transition matrix. Instead this is the matrix $$A$$, since it transfers the state at time $$k$$ to the state at time $$k+1$$. In general, for the discrete linear system $$x_{k+1} = A_kx_{k}$$ you can write the matrix $$\Phi(k,k') = \prod_{j=k'}^{k-1}A_j$$ so that $$x_k = \Phi(k,k')x_{k'}.$$ Note that in the time-invariant case we have $$\Phi(k,k') = A^{k-k'}$$ and $$\Phi(k,0) = A^k$$, which is what you have. Sometimes $$\Phi(k,k')$$ is referred to as the state transition matrix and $$A$$ is referred to as its generator.
The best way to have matlab simulate the solution to $$x_{k+1} = A_kx_k$$, is to use the equation itself rather than trying to calculate any transition matrix (other than for a single step). You can do this by storing the values of the state vector and computing the next value using the update equation. If you really need an analytic formula for $$\Phi(k,k')$$ or $$A^k$$, you can find the Jordan form and take powers of the Jordan blocks, which have well-known formulas that depend on the block type.
• I dont really get your general comments tbh: In the time invariant case, the state transition matrix is clearly given by $\Phi(k, 0)$ (as you wrote), which is equal to $A^k$. How can $A^k$ be "instead [...] the matrix $A$"? Regarding programming, I will check the Jordan form. – SampleTime Mar 29 at 23:45
• @SampleTime $A^k$ is not the same as $A$ in general. You are referring to $A^k$ as the state transition matrix, but this is not usually how the STM is defined. Usually it is either defined as $A$ or as $\Phi(k,k')$. I would look in your text and make sure you are getting the right definition – Mortified Through Math Mar 30 at 1:25
• Ok, I have never seen the use of STM for the $A$ matrix. In every book I have read so far, $A$ is refered to as system matrix while the STM is always defined in terms of the matrix exponential (continous) or power (discrete) of said system matrix. Do you have a reference for that usage? Would be interesting. – SampleTime Mar 30 at 11:22
• @SampleTime I googled "state transition matrix" and got this ocw.mit.edu/courses/electrical-engineering-and-computer-science/… which uses $\Phi$. Yes, sometimes people call $A$ the system matrix, but that is more common in continuous time when it is only an infinitesimal generator of the semigroup. – Mortified Through Math Mar 30 at 12:38