Discretizing the Action for Feynman Path Integral I am trying to understand how to compute path integrals, given a Lagrangian. I understand how it is done for the free particle, but I am confused for other actions. I am having trouble understanding how to decompose the action into subintervals. I understand how to do this for the free particle:
$$S=\frac{m}{2}\frac{(x_b-x_a)^2}{t_b-t_a}\to\frac{m}{2}\sum_{i=1}^N\frac{(x_{i+1}-x_i)^2}{\epsilon}\tag{1}$$
where, $$\begin{align}
\epsilon & = t_{i+1}-t_i \\
N\epsilon & = t_b - t_a \\
t_0  &= t_a \\
t_N & = t_b \\
x_0 &= x_a \\
x_N &= x_b
\end{align}\tag{2}$$
However, I am having trouble understanding how to do this for more general cases. For example, the classical action of the harmonic oscillator is $$S=\frac{m\omega}{2\sin\omega T}\left(\left(x_1^2+x_2^2\right)\cos\omega T -2x_1x_2\right).\tag{3}$$ 
I tried to make this discrete as follows: 
$$\begin{align}S=&\frac{m\omega}{2\sin\omega T}\left(\left(x_1^2+x_2^2\right)\cos\omega T -2x_1x_2\right)\cr
\to&\frac{m\omega}{2\sin\omega\epsilon}\sum_{i=1}^N\left(\left(x_{i+1}^2+x^2_i\right)\cos\omega\epsilon-2x_{i+1}x_i\right).\end{align}\tag{4}$$
Therefore, the path integral becomes 
$$\begin{align}U(x_b,t_b,x_a,t_a)=&\lim_{\epsilon\to 0}\frac{1}{A}\int\cdots\iint\exp\left(\frac{i}{\hbar}\frac{m\omega}{2\sin\omega\epsilon}\sum_{i=1}^N\left(\left(x_{i+1}^2+x^2_i\right)\cos\omega\epsilon-2x_{i+1}x_i\right)\right)\cr
&\times\frac{\mathrm{d}x_1}{A}\frac{\mathrm{d}x_2}{A}\cdots\frac{\mathrm{d}x_{N-1}}{A}\tag{5}\end{align}$$
$$A=\left(\frac{2\pi i\hbar\epsilon}{m}\right)^{\frac{1}{2}}.\tag{6}$$
However, when I evaluated the integral for $x_1$, I got the following:
$$\sqrt{\frac{m\sin\omega\epsilon}{2i\pi\hbar~2\epsilon^2\omega\cos\omega\epsilon}}\exp\left(\frac{i\omega m}{2\hbar  \sin(2\omega\epsilon)}\left(\left(x_0^2+x_2^2\right)\cos2\omega\epsilon-2x_0x_2\right)\right).\tag{7}$$
Comparing to the actual propagator, which is below, I see I got some things right, but the normalization factor is all wrong. Where did I go wrong?
$$U(x_b,t_b,x_a,t_a)=\sqrt{\frac{m\omega}{2\pi i\hbar\sin\omega T}}\exp\left(\frac{im\omega}{2\hbar\sin\omega T}\left(\left(x_a^2+x_b^2\right)\cos\omega T -2x_ax_b\right)\right).\tag{8}$$
P.S. I am using Quantum Mechanics and Path Integrals, Feynman and Hibbs to learn this.
 A: The action that goes in the path integral is a functional of the path $x(t)$. Namely,
$$S[x(t)]=\int_0^T dt \frac{m}{2}\left(\dot{x}^2-\omega^2 x^2\right)$$
which can be discretized very simply
$$S=\frac{m}{2}\sum_{i=1}^N\frac{(x_{i+1}-x_i)^2}{\epsilon}-\epsilon\omega^2 x_i^2$$
The $x$ in this action is not necessarily a classical trajectory, but this is still the action which appears in classical mechanics and leads to the equation of motion of the harmonic oscillator through the Euler-Lagrange equation.
You were confusing this action functional with the numerical value of the action functional evaluated at the classical path
$$S[x_{cl}]=\frac{m\omega}{2\sin\omega T}\left(\left(x_1^2+x_2^2\right)\cos\omega T -2x_1x_2\right)$$
This is just a number (albeit a number that depends on boundary conditions). There is no $x(t)$ appearing in it so you can't vary it to get Euler-Lagrange equations, and you can't integrate it in a path integral.
A: *

*As far as the single Gaussian $x_1$ integration goes, OP is doing the correct thing (up to possible typos) in eq. (7). However since the time-increment $$\epsilon~\ll~ \omega^{-1}\tag{i}$$ is supposed to be small (in order for Feynman's fudge factor $1/A$ to be valid), we have under the square root
$$ \frac{\sin (\omega \epsilon)}{2\epsilon^2\omega\cos (\omega \epsilon)}~\approx~\frac{1}{2\epsilon}~\approx~\frac{\omega }{\sin (\omega 2\epsilon)}, \tag{ii}$$
in agreement with Feynman's formula (8) for the harmonic oscillator. See also this related Phys.SE post.

*It is possible to extend the result from 2 to $N$ time-increments, see e.g. this Phys.SE post. Feynman's formal conceptional idea is to replace the off-shell action $I[x;t_a;t_b]$ in the path integral with the discretized sum of on-shell actions $\sum_{i=1}^N S(x_i,t_i;x_{i-1},t_{i-1})$; then integrate over all intermediate positions $x_1, \ldots, x_{N-1}$; and in the end take the continuum limit $N\to\infty$. 
