Duhamel's Principle for Integrals I am not talking about Duhamel's Integral. That is something different.
There is a lot of coverage of Duhamel's Principle in the context Diff Eq, but I have a copy of "Advanced Calculus" by John M. H. Olmsted copyright 1951. In Chapter 8 page 238 there is a section on Duhamel's Principle for Integrals. It is a weird integral of a function of two variables where two different system of tags of a net are used for the two different arguments to the function.
Duhamel's Principle for Integrals. Let $f\left(t\right)$ and $g\left(t\right)$
be integrable on $\left[a,b\right]$ and let $\phi\left(x,y\right)$
be everywhere continuous. Then, in the sense of section $401$, the limit of the sum $\sum_{i=1}^{n}\phi\left(f\left(t_{i}\right),g\left(t\acute{_{i}}\right)\right)\triangle t_{i}$ where $a_{i-1}\leq t_{i}\leq a_{i}$ and $a_{i-1}\leq t_{i}^{\acute{}}\leq a_{i}$ as the norm of the net  $\Re:a=a_{0},a_{1},\cdots a_{n}=b$ tends toward
zero, exists and is equal to the definite integral  $\intop_{a}^{b}\phi\left(f\left(t\right),g\left(t\right)\right)\,dt$,
which also exits.
\begin{equation}
\lim_{\left|\Re\right|\rightarrow0}\sum_{i=1}^{n}\phi\left(f\left(t_{i}\right),g\left(t\acute{_{i}}\right)\right)\triangle t_{i}=\intop_{a}^{b}\phi\left(f\left(t\right),g\left(t\right)\right)dt
\end{equation}
The point to note in the definition in Olmsted is that $t_{i}$ and $t\acute{_{i}}$ do not have to be equal, just in the same sub-interval.
Olmsted proves Duhamel's principle for the continuous case. Olmstead uses Duhamel's Principle to prove the general formula for the length of an arc.
I have a copy of Baby Rudin, and he proves the same formula for the length of an arc without mentioning Duhamel's Principle for Integrals. When I search the Internet for Duhamel's Principle, all the references to Duhamel's Principle are to the Diff Eq version.
Questions: 
1) Has Duhamel's Principle for Integrals gone out of style?
2) Where can I find a proof of the general case of Duhamel's Principle for Integrals?
3) What is the relationship between Duhamel's Principle for Integrals and the Diff EQ version of Duhamel's Principle?
 A: The proof of Duhamel Principle is based on the following criterion of Riemann integrability given by Riemann himself:

Criterion for Riemann Integrability: Let the function $f:[a, b] \to\mathbb {R} $ be bounded on $[a, b] $. Then $f$ is Riemann integrable on $[a, b] $ if and only if for any given positive real numbers $\epsilon, \sigma$ we can find a positive real number $\delta$ such that for any partition (net) $P$ of $[a, b] $ with norm (mesh) $||P||<\delta$ the subintervals of $P$ on which the oscillation of $f$ is at least $\sigma$ have a total combined length less than $\epsilon$.

The above theorem is a precursor to the famous criterion given by Lebesgue which says that a bounded function is Riemann integrable on a closed interval if and only if the set of its discontinuities on this interval is of measure zero.
Let's assume that functions $f, g$ from $[a, b] $ to $\mathbb {R} $ are Riemann integrable on $[a, b] $ and the function $\phi:\mathbb{R} ^2\to\mathbb {R} $ is continuous on $\mathbb {R} ^2$. Then using the above criterion or Lebesgue's criterion we can show that the function $F:[a, b] \to\mathbb{R} $ defined by $F(x) =\phi(f(x), g(x)) $ is Riemann integrable on $[a, b] $.
Let's start with an arbitrary $\epsilon >0$. Let $A, B$ be positive upper bounds for $|f|, |g|$ on $[a, b] $ respectively and let $M$ be a positive upper bound for $|\phi|$ on the rectanglular region $\mathcal{R} $ with diagonal points $(-A, - B), (A, B) $. We are supposed to find a $\delta>0$ such that for all partitions $P=\{x_0,x_1,x_2,\dots ,x_n\} $ of $[a, b] $ with norm $||P||<\delta $ we have $$\left|\sum_{i=1}^{n}\phi(f(t_i),g(t'_i))(x_i-x_{i-1})-\int_{a} ^{b} \phi(f(x), g(x)) \, dx\right|<\epsilon \tag{1}$$ for any two sets of tags $t_i, t'_i\in[x_{i-1},x_i]$.
To find such a $\delta$ we note that there is a positive $\delta_1$ such that $$\left|\sum_{i=1}^{n}\phi(f(t_i),g(t_i))(x_i-x_{i-1})-\int_{a}^{b}\phi(f(x),g(x))\,dx\right|<\frac{\epsilon} {2}\tag{2}$$ for all partitions $P$ of $[a, b] $ with norm less than $\delta_1$ and any set of tags $t_i\in[x_{i-1},x_i]$. Our job is done if we can find a $\delta_2>0$ such that $$\left|\sum_{i=1}^{n}\{\phi(f(t_i),g(t_i))-\phi(f(t_i),g(t'_i))\}(x_i-x_{i-1})\right|<\frac{\epsilon} {2}\tag{3}$$ for all partitions $P$ with norm less than $\delta_2$ and any sets of tags $t_i, t'_i\in[x_{i-1},x_i]$. Clearly we can choose $\delta=\min(\delta_1,\delta_2)$ and if $P$ is any partition with norm less than $\delta $ then both inequalities $(2),(3)$ hold simultaneously and imply $(1)$.
Thus our desired goal is to ensure the inequality $(3)$ and changing symbols a bit it is sufficient to prove that for any given $\epsilon>0$ there is a $\delta>0$ such that $$\left|\sum_{i=1}^{n}\{\phi(f(t_i),g(t_i))-\phi(f(t_i),g(t'_i))\}(x_i-x_{i-1})\right|<\epsilon\tag{4}$$ for all partitions $P$ with norm less than $\delta$ and any sets of tags $t_i, t'_i\in[x_{i-1},x_i]$.
To find such a $\delta$ we note that $\phi$ is uniformly continuous on the rectanglular region $\mathcal{R} $ and hence there is a $\delta'>0$ such that $$|\phi(x, y) - \phi(x', y') |<\frac{\epsilon} {2(b-a)} $$ for all points $(x, y) $ and $(x', y')$ in $\mathcal{R} $ with $\sqrt{(x-x') ^2+(y-y')^2}<\delta'$. Since $g$ is Riemann integrable on $[a, b] $ there is a positive $\delta$ such that for all partitions $P$ of $[a, b] $ with norm less than $\delta$ the combined length of subintervals of $P$ where the oscillation of $g$ is at least $\delta'$ is less than $\epsilon/4M $. Next consider the sum, henceforth denoted by $S$, on left in the inequality $(4)$ which can be split into two parts based on the subintervals created by partition $P$. The first part, say $S_1$, is based on those sub-intervals $[x_{i-1},x_i]$ where the oscillation of $g$ is less than $\delta'$. For these intervals the difference $$|\phi (f(t_i), g(t_i)) - \phi(f(t_i), g(t'_i)) |<\frac{\epsilon} {2(b-a)}$$ and hence $|S_1|<\epsilon/2$. The second part, say $S_2$, is based on subintervals where the oscillation of $g$ is at least $\delta'$ and the combined length of such subintervals is less than $\epsilon /4M $ so that $|S_2|<\epsilon /2$ and then $$|S|=|S_1+S_2|\leq |S_1|+|S_2|<\epsilon $$
