On page 93 of Mechanics by Landau and Lifshitz, they discuss a general derivation of Kapitza's method. At the end of their derivation, they say that their results can easily be generalised to any number of degrees of freedom, and they quote the following result (see link for context): $$U_{\text{eff}}= U+\frac{1}{2\omega^2}\sum_{ik}a_{ik}^{-1}\overline{f_if_k}.$$ I spent a few hours trying to reproduce the result, but always got stuck. This is what I have tried so far:
Let's assume the kinetic energy of the system has the form: \begin{equation} \label{eq:kinener} T = \frac{1}{2}\sum_k \sum_j t_{jk}\dot{q_j}(t)\dot{q_k}(t) \end{equation} where the sum is over all the degrees of freedom. The Euler-Lagrange equations lead to \begin{equation} \label{bwv} \sum_k t_{jk}\ddot{q_k}(t) = -\frac{\partial U(\mathbf{q}(t))}{\partial q_j} + f^{(j)}(\mathbf{q}(t)), \,\,\,\,\,\,\, 0 \leq j \leq n, \end{equation} where I added a "generalised force disturbance" $f^{(j)}$ . Let's now separate the motion into fast and slow components: \begin{equation} \mathbf{q}(t) = \mathbf{Q}(t) + \boldsymbol{\xi}(t), \end{equation} where $\boldsymbol{\xi}(t)$ corresponds to the small oscillations.
Using this, and expanding to first order around the smooth component, I get:
\begin{multline} \label{eq:bwvt} \sum_k t_{jk}(\ddot{Q_k}+ \ddot{\xi_k}) = -\frac{\partial U(\mathbf{Q})}{\partial q_j} - \sum_k \frac{\partial^2 U(\mathbf{Q})}{\partial q_k\partial q_j} \xi_k \\ + f^{(j)}(\mathbf{Q},t) + \sum_k\xi_k\frac{\partial f^{(j)}(\mathbf{Q})}{\partial q_k}, \end{multline}
We now fix $\boldsymbol{\xi}(t)$ by \begin{equation} \label{eq:choiceXi} \sum_k t_{jk}\ddot{\xi_k}(t) = f^{(j)}(\mathbf{Q},t). \end{equation} Integrating, I regard $\mathbf{Q}$ constant and discard all integration constants: \begin{equation} \sum_k t_{jk} \xi_k(t) = -\frac{1}{\omega^2}f^{(j)}(\mathbf{Q},t) \implies \xi_k = -\frac{1}{\omega^2}\sum_j t_{kj}^{-1}f^{(j)}(\mathbf{Q},t), \end{equation} where $t_{kj}^{-1}$ are the inverse matrix elements. Filling this all in and averaging: \begin{align*} \sum_k t_{jk} \ddot{Q}_k &= -\frac{\partial U(\mathbf{Q})}{\partial Q_j} + \sum_k \langle \xi_k \frac{\partial f^{(j)}(\mathbf{Q},t)}{\partial Q_k} \rangle \\ &= -\frac{\partial U(\mathbf{Q})}{\partial Q_j} -\frac{1}{\omega^2} \sum_k\sum_i t_{ki}^{-1} \langle f^{(j)}(\mathbf{Q},t) \frac{\partial f^{(j)}(\mathbf{Q},t)}{\partial Q_k} \rangle \\ &=-\frac{\partial U(\mathbf{Q})}{\partial Q_j} -\frac{1}{2\omega^2} \sum_k\sum_i t_{ki}^{-1} \langle \frac{\partial f^{(j)2}(\mathbf{Q},t)}{\partial Q_k}\rangle \end{align*} ... this is where I'm stuck...
Could it be that the Landau and Lifshitz result is just wrong?