A mathematical motivation for Lax-Friedrich type of Numerical Fluxes A Lax-Friedrichs (LF) type of flux for a conservation law $\partial_tU+\partial_xf(U)=0$ is given by 
\begin{align}
F(U^-, U^+) = \frac{1}{2} \Big(f(U^-) + f(U^+)\Big)\cdot \nu - \frac{1}{2} \lambda(U^+ - U^-)
\end{align}
where for a classical LF flux in one dimension, $\lambda  = \Delta t / \Delta x$. I understand that the motivation comes from trying to write the central flux scheme 
\begin{align}
U^{n+1} = U^{n} - \frac{\Delta t}{\Delta x} 
\Big(f(U^n_{i+1}) - f(U^n_{i-1})\Big)
\end{align}
by 
\begin{align}
U^{n+1} = \frac{1}{2}\left(U^n_{i+1} + U^n_{i-1} \right)
- \frac{\Delta t}{\Delta x} \Big(f(U^n_{i+1}) - f(U^n_{i-1})\Big)
\end{align}
which then covers a wider domain than before and hence is more stable.
The step of replacing the value at a point by average seems like a small trick. I am curious what is the idea behind such scheme, and if we can arrive at the scheme by writing a modified version of PDE and proceeding along more natural path.
 A: The most natural way to derive Lax-Friedrichs's scheme is to consider the discretisation where the approximate state $U$ is constant ($\equiv U_i^n$) in cells $((i-1)\Delta x,(i+1)\Delta x)\times((n-1)\Delta t,n\Delta t)$ where $n+i\in2{\mathbb Z}$. To pass from the array $(U_i^n)_i$ to the next one $(U_i^{n+1})_i$, you start from the piecewise constant data at time $n\Delta t$~;  you can solve the Cauchy problem explicitly by gluing the solutions of the Riemann problems (a RP is a Cauchy problem with data constant on each side of one point). These RP do not interact on a time interval $\Delta t$ guaranted by the Courant-Friedrichs-Levy inequality. The solution at $(n+1)\Delta t$ is projected over the space of piecewise constant states (consisting in taking mesh averages).
If instead you consider all the cells $((i-\frac12)\Delta x,(i+\frac12)\Delta x)\times((n-1)\Delta t,n\Delta t)$ ($i,n\in{\mathbb Z}$), the same procedure yields the Godunov's flux, which is less explicit, though still computable.
The stability of these procedures is quite obvious, because both stages (explicit solution on a time step, projection) is stable.
