Non-constant coefficient matrix in first order linear differential equations I want to solve a differential equation of the following form 
$$
\frac{d}{dt}x=A(t)x\, ,
$$
where $A(t)$ does not commute at different times. This equation holds on the interval $(a,b)$. Hence, the solution cannot be written as $e^{\int_a^td\tau A(\tau)}$. I have been searching now for several hours, and I only found (source: http://www.macs.hw.ac.uk/~simonm/linalg.pdf pg. 92) the so called "Neumann series" 
$$
x(t)=\left(I+\int_0^{t}A(\tau)d\tau + \int_0^{t}A(\tau_1)\int_0^{\tau_1}A(\tau_2)d\tau_2 d\tau_1+...\right) x_0\, .
$$
(The series is very reminiscent of the Dyson series with the evolution operator in quantum mechanics, here $A(t)$ is not hermitian.)
Here are my doubts:


*

*In the formula the lower bound is $0$. Is this part of the formula or can it be any number? in my case $a$.

*Is this an infinite sum? if yes how can I stop at a certain summand to have an approximation?

*How does this series continue? is it similar to the time ordering in
quantum mechanics, i.e., one takes every possible commutation or is the third term just: 
$$
 \int_0^{t}A(\tau_1)\int_0^{\tau_1}A(\tau_2)\int_0^{\tau_2}A(\tau_3)d\tau_2d\tau_3 d\tau_1\, .
$$

*In the link above, it says, one can check that this is the solution by taking the derivative with respect to $t$. How can one take the derivative, when it depends on the upper bound of the integral and how is the product rule for non commutative matrices?


Thank you very much.
 A: I think you probably already found the answers to your questions, but I'll still provide one for the others who like you and me would be interested by that:

*

*Like for any Cauchy problem, $\,\,\dot{\vec{y}}(t)=f(t,\vec{y}(t))\,\,$, you must specify an initial value $(t_0,\vec{y}_0)$ that the solution $\vec{y}(t)$ should satisfy to. In other words, this should hold true: $\vec{y}(t_0)=\vec{y}_0$. 
So as you might see, the author chooses the initial value to be $(0,\vec{y}_0)$ and with the Ansatz that their propose this makes totally sense since:
$$ \vec{y}(0)=\left(I+\int_0^0A(\tau)\,d\tau+\int_0^0 A(\tau_1)\int_0^{\tau_1}A(\tau_2)\,d\tau_2\,d\tau_1\,+\,...\right)\vec{y}_0 = I\,\vec{y}_0=\vec{y}_0$$
So if now your initial value is $(t_0,\vec{y}_0)$, then the Ansatz for your solution would be:
$$ \vec{y}(t)=\left(I+\int_{t_0}^t A(\tau)\,d\tau+\int_{t_0}^t A(\tau_1)\int_{t_0}^{\tau_1}A(\tau_2)\,d\tau_2\,d\tau_1\,+\,...\right)\vec{y}_0$$
such that when you substitute $t=t_0$, the integrals vanishe and you get your desired equality $\vec{y}(t_0)=\vec{y}_0$. 
Replace $t_0$ by $a$ and you will obtain the formula you want.


*As the author mentions it, this Ansatz for the solution of the initial value problem is written as a Neumann series which is a generalization of geometric series to operators, here the operator is:
$$
T[O]=\int_{t_0}^{t} A(\tau)\,O(\tau)\,d\tau
$$
with $O$ being a matrix of the same size as A. Then the Ansatz can be rewritten as:
$$
\vec{y}(t)= \sum_{k=0}^\infty T^k[\,I\,] \,\, \vec{y}_0 \qquad \mathrm{where}
\qquad T^k[\,I\,] := T^{k-1} \circ T[\,I\,] \equiv T^{k-1}[\,T[\,I\,]\,]
$$
and $T^0[\,I\,] \equiv I\,$. This particular geometric series has also a convergence criteria which is that $|T|<1$ in the operator norm. We must typically prove that $T$ is a bounded linear operator on a Banach space. Then you will probably see that $|T|<1$ is equivalent to $\int_{t_0}^t ||A(\tau)||\,d\tau\,< \infty\,$. $\quad$ [to be done] 
So to answer your question: If your matrix $A$ satisfies the convergence criteria, then yes stopping at a particular term in the sum above will give a more or less good approximation of your solution. Typically higher terms will contribute less than the first ones.


*Yes, this is the third term of the series. Just take the operator $T$ defined above and repeatedly apply it to the last term to obtain the next one.


*As the author mentions this is straightforward to prove. I use in the following lines the Leibniz integral rule since a matrix integration is simply a component wise integration.
$$
\begin{align}
\frac{d}{dt}\vec{y}(t) &= \frac{d}{dt}\left[ \left(I+\int_{t_0}^t A(\tau)\,d\tau+\int_{t_0}^t A(\tau_1)\int_{t_0}^{\tau_1}A(\tau_2)\,d\tau_2\,d\tau_1\,+\,...\right)\vec{y}_0 \right] = \\
&= \left(0+\frac{d}{dt}\left[\int_{t_0}^t A(\tau)\,d\tau \right] + \frac{d}{dt}\left[\int_{t_0}^t A(\tau_1)\int_{t_0}^{\tau_1}A(\tau_2)\,d\tau_2\,d\tau_1\,\right]+\,...\right)\vec{y}_0 = \\
&=\left(A(t)\,\frac{d}{dt} \bigg[\,t\,\bigg] + A(t)\,\int_{t_0}^t A(\tau_2)\,d\tau_2 \,\frac{d}{dt} \bigg[\,t\,\bigg] + \,...\right)\vec{y}_0 = \\
&= A(t) \,\, \left(I+\int_{t_0}^t A(\tau)\,d\tau+\int_{t_0}^t A(\tau_1)\int_{t_0}^{\tau_1}A(\tau_2)\,d\tau_2\,d\tau_1\,+\,...\right)\vec{y}_0 =\\
&= A(t)\,\vec{y}(t)
\end{align}
$$
So with that we have proven that the Ansatz is indeed a solution to the differential ordinary matrix equation. However, I must emphasize once again, the problem must satisfy the convergence criteria of the Neumann series. 
(I did not get the question about the product rule, sorry)
