∏

How to Find the Product Integral?

Notational Conventions

In the subsequent text we use the following convention: A function $A:[a,b] \to\real^{n\times n}$ is called Riemann integrable, if its entries $a_{ij}$ are Riemann integrable functions on $[a,b]$. In this case we put

$$ \int_a^b A(t)\,\d t=\left\{\int_a^b a_{ij}(t)\,\d t\right\}_{i,j=1}^n . $$

We will often encounter integrals of the type

$$ \int_a^b\!\!\int_a^{x_k}\!\!\cdots\!\!\int_a^{x_2} A(x_k)\cdots A(x_1)\,\d x_1\cdots\,\d x_k. $$

These integrals should be interpreted as iterated integrals, where $x_k\in[a,b]$ and $x_i\in[a,x_{i+1}]$ for $i\in\{1,\ldots,k-1\}$.

Expansion into Peano Series

Theorem.

If $A:[a,b]\to\real^{n\times n}$ is a Riemann integrable function, then both product integrals exist and

$$ \prod_a^b (I+A(x)\,\d x)=I+\sum_{k=1}^\infty\;\int_a^b\!\!\int_a^{x_k}\!\!\cdots\!\!\int_a^{x_2} A(x_k)\cdots A(x_1)\,\d x_1\cdots\d x_k, $$
$$ (I+A(x)\,\d x)\prod_a^b =I+\sum_{k=1}^\infty\;\int_a^b\!\!\int_a^{x_k}\!\!\cdots\!\!\int_a^{x_2} A(x_1)\cdots A(x_k)\,\d x_1\cdots\d x_k. $$
Proof:

Volterra's proof goes as follows: Expanding the product $P(A,D)$ gives

$$ \prod_{i=m}^1(I+A(\xi_i)\Delta t_i)= I+\sum_{k=1}^m\left(\sum_{1\leq i_1<\cdots<i_k\leq m}\!\! A(\xi_{i_k})\cdots A(\xi_{i_1})\Delta t_{i_1}\cdots \Delta t_{i_k}\right). $$
Where does this complicated sum come from? (Expanding the product)

Don't let the multi-index notation scare you! It arises from simple polynomial expansion. Imagine you are multiplying terms formatted as $(I + A_i)$.

If we expand a product of just two brackets, we get: $$ (I + A_2)(I + A_1) = I + A_1 + A_2 + A_2 A_1 $$

If we add a third one, the product starts building up and gets more complex: $$ (I + A_3)(I + A_2)(I + A_1) = I + (A_1 + A_2 + A_3) + (A_2 A_1 + A_3 A_1 + A_3 A_2) + A_3 A_2 A_1 $$

Notice the emerging pattern: We always get the identity matrix $I$, then a sum of single elements, then a sum of all possible pairs (ordered from highest index to lowest), then triples, and so on.

The condition $1 \leq i_1 < \cdots < i_k \leq m$ in the theorem is just a rigorous mathematical way of saying: "Choose all possible combinations of $k$ matrices, keeping their relative order from the original product."

Volterra now argues that for $\nu(D)\to 0$ we obtain

$$ \sum_{1\leq i_1\leq m}A(\xi_{i_1})\Delta t_{i_1}\to $$
$$ \int_a^b A(x_1)\,\d x_1, $$
$$ \sum_{1\leq i_1<i_2\leq m}A(\xi_{i_2})A(\xi_{i_1})\Delta t_{i_1}\Delta t_{i_2}\to $$
$$ \int_a^b\!\!\int_a^{x_2} A(x_2)A(x_1)\,\d x_1\,\d x_2, $$

and generally

$$ \sum_{1\leq i_1<\cdots<i_k\leq m}\!\!\!\!\!\!\!\! A(\xi_{i_k})\cdots A(\xi_{i_1})\Delta t_{i_1}\cdots \Delta t_{i_k}\to $$
$$ \int_a^b\!\!\int_a^{x_k}\!\!\cdots\!\!\int_a^{x_2} A(x_k)\cdots A(x_1)\,\d x_1\cdots\d x_k $$

for every $k\in\{1,\ldots,m\}$. Using the fact that $m\to\infty$ for $\nu(D)\to0$, Volterra arrived at the result

$$ \prod_a^b (I+A(x)\,\d x)=I+\sum_{k=1}^\infty\;\int_a^b\!\!\int_a^{x_k}\!\!\cdots\!\!\int_a^{x_2} A(x_k)\cdots A(x_1)\,\d x_1\cdots\d x_k. $$

Special Cases & Solvable Forms

1. Scalar Functions

Example:

If $n=1$, then the function $A:[a,b]\to\real^{n\times n}$ is in fact a scalar function, and we usually write $\prod_a^b(1+A(t)\,\d t)$ instead of $\prod_a^b(I+A(t)\,\d t)$.

Using Theorem [volt-peano] and Lemma [volt-commute-lemma] we obtain

$$ y(x)=\prod_a^x(1+A(t)\,\d t)=1+\sum_{k=1}^\infty\frac{1}{k!}\left(\int_a^x A(t)\,\d t\right)^k=\exp\left(\int_a^x A(t)\,\d t\right), $$

which is indeed a solution of the differential equation $y'(x)=A(x)y(x)$ and satisfies $y(a)=1$.

Scalar case solved using logarithmic differentiation

An alternative and very intuitive way to solve the scalar equation $y'(x) = A(x)y(x)$ with $y(a)=1$ is via standard calculus. Assuming $y(x) > 0$, we can divide by $y(x) $ to get: $$ \frac{y'(x)}{y(x)} = A(x) $$

Recognizing the left side as the derivative of the natural logarithm gives $\frac{\mathrm{d}}{\mathrm{d}x}(\ln y(x)) = A(x)$. Integrating both sides from $a$ to $x$ yields: $$ \ln y(x) - \ln y(a) = \int_a^x A(t)\,\mathrm{d}t $$

Since $y(a) = 1$ and $\ln(1) = 0$, taking the exponential of both sides immediately leads to the solution $y(x) = \exp\left(\int_a^x A(t)\,\mathrm{d}t\right)$.

Crucial Warning: This elegant logarithmic method is strictly limited to the scalar case ($n=1$). If $A(x) $ and $Y(x)$ are matrices, the derivative of a matrix logarithm does not obey the simple scalar rule. Specifically: $$ \frac{\mathrm{d}}{\mathrm{d}x} \Big(\ln Y(x)\Big) \neq Y'(x)Y(x)^{-1} $$

...unless the matrix $Y(x)$ commutes with its own derivative ($Y(x)Y'(x) = Y'(x)Y(x)$), which is rarely true for time-dependent matrix systems. This is exactly why Volterra had to invent the Product Integral and the Peano series expansion in the first place—to handle the non-commutative matrix cases where logarithms fail us!

2. Constant Matrix Functions

Example:

Recall that if $A\in\real^{n\times n}$, then the exponential of $A$ is defined as

$$ \exp A=\sum_{k=0}^\infty\frac{A^k}{k!}. $$

The fundamental matrix of the system of equations

$$ y_i'(x)=\sum_{j=1}^n a_{ij}y_j(x),\;\;\; i=1,\ldots,n $$

is given by

$$ Y(x)=\prod_a^x(I+A\,\d t)=I+\sum_{k=1}^\infty\;\int_a^b\!\!\int_a^{x_k}\!\!\cdots\!\!\int_a^{x_2} A^k\,\d x_1\cdots\d x_k= $$
$$ I+\sum_{k=1}^\infty\frac{(x-a)^kA^k}{k!}=e^{(x-a)A} $$

We also remark that a similar calculation leads to the relation

$$ (I+A\,\d t)\prod_a^x=e^{(x-a)A},\;\;\; x\in[a,b]. $$