Free energy of coupled classical harmonic oscillators I'm looking to find the thermodynamic (NVT) free energy of a classical coupled harmonic oscillator system such as the one below:

(image taken from http://openmetric.org/StatisticalPhysics/equilibrium/week3.html)
I would like a solution that allows an arbitrary $N$ number of masses, and ideally I would like to have a fully general expression with arbitrary (different) spring constants and arbitrary (different) masses.
I tried to compute this free energy by hand by computing the partition function:
$$ Z = \int_{-\infty}^{\infty}d\vec{p} \int_{-\infty}^{\infty}d\vec{x}\; e^{-\beta H}$$
where 
$$ H = \sum_{i=1}^{N} \frac{p_i^2}{2m_i} + \sum_{i=0}^{N} m_i\omega_i^2(x_{i+1}-x_i)^2 $$
and $x_i$ denotes the displacement from equilibrium of the $i$th block, and $x_0=x_{N+1}=0$ represent the walls at the ends.
I was able to derive the expressions for the free energy $F = -\frac{1}{\beta}\ln Z$ for one, two, and three blocks with identical masses and identical springs (with the hope of seeing an extendible pattern) but sadly no obvious patterns emerged. The calculations are also rather tedious.
I don't doubt that this has already been done many times before -- does anyone have a reference to a solution?
 A: Make the change of variable to $\delta_i = x_{i+1} - x_i,\;  i=0\dots N-1.$ Then the system is uncoupled and $Z = \prod_i z_i$ with
$$ z_i = \int e^{-\beta p^2/2m_i}dp \times \int e^{-\beta m_i\omega_i^2\delta^2}d\delta = \frac{\pi\sqrt{2}}{\beta\omega_i}.$$
A: Let $k_i = m_i \omega_i^2$, and we have
\begin{equation}
H\big(\{p_i\},\{x_i\}\big) = \sum_{i=1}^{N} \frac{p_i^2}{2m_i} + \sum_{i=0}^{N} \frac{k_i}{2}(x_{i+1}-x_i)^2,
\end{equation}
where $x_{{\scriptscriptstyle N+1}} = x_{{\scriptscriptstyle 0}} = 0$. The partition function is given by
\begin{equation}
\begin{split}
Z &= \int \left(\prod_{i=1}^N dp_i\right) \left(\prod_{i=1}^N dx_i\right) \,\exp\left[-\frac{H\big(\{p_i\},\{x_i\}\big)}{k_{B}T}\right]\\
&= \int \left(\prod_{i=1}^N dp_i\right)\, \exp\left(-\frac{1}{k_{B}T}\sum_{i=1}^{N} \frac{p_i^2}{2m_i}\right)\int \left(\prod_{i=1}^N dx_i\right)\, \exp\left[-\frac{1}{k_{B}T}\sum_{i=0}^{N}\frac{k_{i}}{2}(x_{i+1} - x_{i})^2\right]\\
&=I_p I_x\,,
\end{split}
\end{equation}
where $I_p$ and $I_x$ respectively represent the $p$ and $x$ integrals in the second line.
The $p$ integral is a product of Gaussian integrals, which gives
\begin{equation}
I_p = (2\pi k_B T)^{N/2}\sqrt{\prod_{i=1}^{N}m_{i}}
\end{equation}
To evaluate $I_x$, we first make the change of variables $x_{i}^\prime = x_{i} / \sqrt{k_B T} $. Then, $I_x = (k_B T)^{N/2} \tilde{I}_{x}$, where
\begin{equation}
\begin{split}
\tilde{I}_{x} &= \int\left(\prod_{i=1}^N dx_i^\prime\right) \,\exp\left[-\sum_{i=0}^{N}\frac{k_{i}}{2}(x_{i+1}^\prime - x_{i}^\prime)^2\right]\\
&= \int\left(\prod_{i=1}^{N+1} dx_i^\prime\right) \,\delta(x_{N+1}^\prime)\,\exp\left[ -\sum_{i=0}^{N}\frac{k_{i}}{2}(x_{i+1}^\prime - x_{i}^\prime)^2\right].
\end{split}
\end{equation}
In fact, for all practical purposes, it suffices to stop at this point. That is,
\begin{equation}
Z = I_p I_x = (2\pi)^{N/2} (k_B T)^N \tilde{I}_x,
\end{equation}
where $\tilde{I}_x$ is independent of any thermodynamic variable of the system (in this case $T$).
Still, it is possible to evaluate the integral $\tilde{I}_x$. To do this, we make another change of variables $y_{i} = x_{i}^\prime - x_{i - 1}^\prime \ \ (i = 1, 2, \ldots, N)$. The associated Jacobian determinant is equal to 1, and it follows that $x_{N+1}^\prime = \sum_{i=1}^{N+1} y_{i}$. Also, we can represent the delta function as
\begin{equation}
\delta(x_{N+1}^\prime) = \int \frac{dq}{2\pi}\, \exp(iqx_{N+1}^\prime)
\end{equation}
Then,
\begin{equation}
\begin{split}
\tilde{I}_x &= \int\frac{dq}{2\pi}\left(\prod_{i=0}^{N} dy_{i+1}\right) \,\exp\big[iq({\scriptstyle \sum_{i=0}^{N} y_{i+1}})\big] \exp\left( -\sum_{i=0}^{N}\frac{k_{i}}{2}y_{i+1}^2\right)\\
&= \int \frac{dq}{2\pi} \prod_{i=0}^{N}\Bigg[\int dy_{i+1} \exp\bigg(-\frac{k_{i}}{2}y_{i+1}^2 + iqy_{i+1}\bigg)\Bigg]\\
&=  \int \frac{dq}{2\pi} \prod_{i=0}^{N} \sqrt{\frac{2\pi}{k_i}}\exp\left(-\frac{q^2}{2k_i}\right)\\
&=\frac{(2\pi)^{N/2}}{\sqrt{\Big(\sum_{i=0}^{N}\frac{1}{k_{i}}\Big)\prod_{i=0}^{N} k_i}}.
\end{split} 
\end{equation}
Collecting the results we have derived so far, we obtain
\begin{equation}
Z = I_p I_x = (2\pi k_B T)^{N} \sqrt{\frac{\prod_{i=1}^{N}m_{i}}{\Big(\sum_{i=0}^{N}\frac{1}{k_{i}}\Big)\prod_{i=0}^{N} k_i}}.
\end{equation}
