Calculation of variance of complicated random variable ( white noise discretization ) so I have been doing some state estimation and in one part of my work it is necessary to discretize a continous time differential equation with white noise. I understood the discretization process for differntial equation ( deterministic part ) but I do not understand the last equation in deriving the covariance of discretized white noise.
Assuming that  A is some matrix and that $x(t)$ is some differentiable vector-valued function mapping into $\mathbb{R}$ function. Further assuming that $w(t,\omega)$ is some continous time wide-sense stationary stohastic process such that $w:\mathbb{R}\times\Omega \to \mathbb{R}^4$ has mean value zero and autocorrelation of all random variables indexed by this process is delta function ( basically said $w(t,\omega)$ is considered white noise ).Further random variables indexed by this process have covaraince matrix Q.Let now :
$$\dot{x}(t) = Ax(t) + w(t)$$
Now using the procedure as described in Derivation of discretization on wikipedia I obtain the following equation:
$$x[k+1]=e^{AT_s}x[k] + \int_{kT_s}^{(k+1)T_s}e^{A[(k+1)T_s -l]}w(l)dl$$
$$u[k]=\int_{kT_s}^{(k+1)T_s}e^{A[(k+1)T_s -l]}w(l)dl$$
where $T_s$ is some discretization constant and $$x[k]=x(kT_s)$$
Now covariances of random variables indexed by $u[k]$ should be :
$$\mathbb{E}[u[k]*u^{T}[k]]=\mathbb{E}[\int_{kT_s}^{(k+1)T_s}e^{A[(k+1)T_s -l]}w(l)dl * \int_{kT_s}^{(k+1)T_s}w^{T}(t)e^{A^{T}[(k+1)T_s -t]}dt]$$
And I am unable to calculate any further. Wherever I look they seem to use the following identity which I would like to be able to justify but cant:
$$\mathbb{E}[\int_{kT_s}^{(k+1)T_s}e^{A[(k+1)T_s -l]}w(l)dl * \int_{kT_s}^{(k+1)T_s}w^{T}(t)e^{A^{T}[(k+1)T_s -t]}dt]=\mathbb{E}[\int_{kT_s}^{(k+1)T_s}\int_{kT_s}^{(k+1)T_s}e^{A[(k+1)T_s -l]}w(l)w^{T}(t)e^{A^{T}[(k+1)T_s -t]}dldt]$$
Could someone justify this for me,or are engineers playing sneaky mathematics again?
EDIT: I understand that in current presentation this is not rigorous and that the differential equation above mathematically has no meaning.Sadly I am not knowledgable in area of Ito Intgration and Stohastic differential equations,so answers avoiding reference to such would be more useful. Material references are welcome
 A: At this point you have the product of two definite integrals that are constants to each other, so you can calculate
\begin{align}
\int_a^b u(s)\,ds\cdot \int_a^bv(t)\,dt
&=\int_a^b u(s)\cdot \left[\int_a^bv(t)\,dt\right]\,ds
\end{align}
Now $u(s)$ is a constant relative to the integral of $v$, so that one can insert this constant into the inner integral
\begin{align}
...&=\int_a^b\left[\int_a^b u(s)\cdot v(t)\,dt\right]\,ds
\end{align}
The interesting part is in the next step where you exchange the integral over the probability space of the expectation value with the two time integrations. The probability argument is usually ommited, putting it back in gives
\begin{align}
{\Bbb E}\left[\int_a^b\int_a^b u(s)\cdot v(t)\,dt\,ds\right]
&=\int_Ω\left[\int_a^b\int_a^b u(\omega,s)\cdot v(ω,t)\,dt\,ds\right]dP(ω)
\end{align}
There you need some Fubini-type result, using compactness of the time intervals and continuity in time of the integrands to obtain boundedness for each $ω\inΩ$, and the finiteness and positivity of the probability measure to get the absolute boundedness of the joint integration over $[a,b]^2\times Ω$. Then
\begin{align}
... 
&=\int_a^b\int_a^b \int_Ω\left[u(\omega,s)\cdot v(ω,t)\right]dP(ω)\,dt\,ds
=\int_a^b\int_a^b {\Bbb E}\left[u(s)\cdot v(t)\right]\,dt\,ds.
\end{align}
