One very useful application of the Cayley-Hamilton theorem is in finding explicit formulas for matrix functions. While this may seem very arbitrary, matrix functions are extremely helpful in solving systems of differential equations.

Consider an n x n square matrix A which has characteristic equation: (-1)^{n}λ^{n}+c_{n-1}λ^{n-1} + ... + c_{1}λ + c_{0} = 0. By the Cayley-Hamilton theorem, the matrix A satisfies it's own characteristic equation so:

(-1)^{n}A^{n}+c_{n-1}A^{n-1} + ... + c_{1}A + c_{0}I = 0.

Which means that A^{n} can be expressed as an (n-1)th degree polynomial function of A and similarly, λ^{n} can be expressed as the same (n-1)th degree polynomial as a function of λ instead of A. So we can find an explicit formula for A^{n} by solving the system of equations which arise for various eigenvalues. In the case of k repeated eigenvalues, the equation for λ^{n} can be differentiated k times to ensure a total of n linearly independent equations so the system may be solved.

We can extend this theory to obtain explicit formulas for other matrix functions like e^{A}, sinA, logA, or even A^{k}.

Consider a function f(x) with Taylor series ∑a_{k}x^{k} which is convergent for all x. By the Cayley-Hamilton theorem, f(λ)=∑a_{k}λ^{k} and f(A) = ∑a_{k}A^{k}. Here the summation goes from 0 to infinity.

Now given that f(A) = ∑a_{n}A^{n}, let's define a function q(A) such that ∑a_{k}A^{k} = q(A)*{(-1)^{n}A^{n}+c_{n-1}A^{n-1} + ... + c_{1}A + c_{0}I}. Where summation on the left goes from n to infinity. Notice that the summation on the right is actually the characteristic equation of the matrix A which is identically zero. Hence we can derive the formula for the function of a matrix f(A) = s_{n-1}A^{n-1} + ... + s_{1}A + s_{0}I similarly f(λ) = s_{n-1}λ^{n-1} + ... + s_{1}λ + s_{0}. For the n x n matrix an explicit formula for the function f(A) can be found by the solving the system of equations for f(λ).

Having such formulas help reduce solving systems of differential equations to something similar to solving one. For example, the system dX/dt = AX can be solved by a number of methods but by comparing it to a linear homogenous first order differential equation, we can immediately notice the solution is the matrix exponential X = e^{At}.

As with matrices, matrix functions may or may not be commutative. In other words, e^{A}e^{B} may not equal e^{A+B} or e^{B}e^{A}. Note also that for functions that are not convergent for all x, all |λ|s must be within the radius of convergence.

In sum, given a matrix and its eigenvalues, it is possible to derive a formula for function of that matrix. The process is very simple and involves solving one system of linear equations. For clarity, let us use a 3 x 3 matrix, A, as an example with eigenvalues λ_{1}, λ_{2}, λ_{3}. Now to find the expression for a function of A, we must first ensure that the eigenvalues are within the radius of convergence of that function. Meaning that the absolute value of the greatest eigenvalue must be less than the radius of convergence of the function's taylor series, otherwise the system may not have solutions. Also, for obvious reasons, f(λ_{n}) must be defined. The functions can be derived by determining the coefficients of the system:

f(λ_{1}) = a_{2}λ_{1}^{2} + a_{1}λ_{1} + a_{0}

f(λ_{2}) = a_{2}λ_{2}^{2} + a_{1}λ_{2} + a_{0}

f(λ_{3}) = a_{2}λ_{3}^{2} + a_{1}λ_{3} + a_{0}

In the case of repeated eigenvalues, we can differentiate the equation to obtain new, distinct equation. Supposing λ_{1}=λ_{2} our system would be:

f(λ_{1}) = a_{2}λ_{1}^{2} + a_{1}λ_{1} + a_{0}

(d/dt)f(λ_{1}) = 2a_{2}λ_{1} + a_{1}

f(λ_{3}) = a_{2}λ_{3}^{2} + a_{1}λ_{3} + a_{0}

And in the case of three repeated eigenvalues:

f(λ_{1}) = a_{2}λ_{1}^{2} + a_{1}λ_{1} + a_{0}

(d/dt)f(λ_{1}) = 2a_{2}λ_{1} + a_{1}

(d^{2}/dt^{2})f(λ_{1}) = 2a_{2}

Solving for a_{2}, a_{1}, and a_{0} we obtain the expression:

f(A) = a_{2}A^{2} + a_{1}A = a_{0}I

Source: Zill and Cullen, Advanced Engineering Mathematics 3rd ed. Jones and Bartlett