The Problem
I currently have the following equation, which I would like to solve using harmonic balancing.
$$\left[\mathbf M_c + \mathbf M_v(\mathbf \Theta)\right]\ddot{\mathbf \Theta} + \mathbf C\dot{\mathbf \Theta} + \left[\mathbf K + \mathbf G(\mathbf \Theta) \right]\mathbf \Theta = \mathbf F(\mathbf \Theta) \tag{1}$$
where $\mathbf M_c$, $\mathbf M_v$, and $\mathbf G$ are all diagonal matrices, and the elements of the latter two are sinusoidal functions.
However, when I compare the results of the harmonic balancing solution with the typical time-domain solution, they are significantly different.
I can get the two solutions to match when I have an equation similar to Eq. 1 above, but without the $\mathbf M_v$ and $\mathbf G$ terms, so that shows me that I have at least some part of the process right. Thus things go wrong when I introduce the variable coefficients $\mathbf M_v$ and $\mathbf G$.
My Attempt
Let us first define the Fourier series in exponential form for the following terms:
$$\mathbf \Theta = \sum^{+\infty}_{n=-\infty}\theta_n e^{in\omega t}, \quad \mathbf F = \sum^{+\infty}_{n=-\infty}\mathbf f_n e^{in\omega t}, \quad \mathbf M_v = \sum^{+\infty}_{k=-\infty}\mathbf m_ke^{ik\omega t}, \quad \mathbf G=\sum^{+\infty}_{k=-\infty}\mathbf g_ke^{ik\omega t}$$
Substituting the above into Eq. 1 I get:
\begin{equation} \left[\mathbf M_c + \sum_k\mathbf m_ke^{ik\omega t} \right]\sum_n (-n\omega)^2\theta_n e^{in\omega t} + \mathbf C\sum_n(in\omega)\theta_ne^{in\omega t} \\ +\left[\mathbf K + \sum_k\mathbf g_ke^{ik\omega t}\right]\sum_n\theta_ne^{in\omega t} = \sum_n\mathbf f_n e^{in\omega t} \tag{2} \end{equation}
Simplifying the above I get:
\begin{equation} \sum_n\sum_k \left[\mathbf g_k -(n\omega)^2\mathbf m_k\right]\theta_ne^{i(k+n)\omega t} + \\ \sum_n\left[in\omega\mathbf C - (n\omega)^2\mathbf M_c + \mathbf K \right]\theta_ne^{in\omega t} = \sum_n\mathbf f_n e^{in\omega t} \tag{3} \end{equation}
So by considering the positive indices of the sums, I turn the above into a matrix of linear equations, by trying to find a common value for the exponential function of all terms, such that it cancels out: e.g. if the exponential with $n$ by itself is equal to 3, then the exponential with $n+k$ means that $n=1,2$ when $k=2,1$:
\begin{align} &\left[i\omega\mathbf C - (\omega)^2\mathbf M_v + \mathbf K\right]\theta_1 = \mathbf f_1\\ &\left[2i\omega\mathbf C - (2\omega)^2\mathbf M_v + \mathbf K\right]\theta_2 + \left[\mathbf g_1 - (\omega)^2\mathbf m_1\right]\theta_1 = \mathbf f_2 \\ &\left[3i\omega\mathbf C - (3\omega)^2\mathbf M_v + \mathbf K\right]\theta_3 + \left[\mathbf g_2 - (\omega)^2\mathbf m_2\right]\theta_1 + \left[\mathbf g_1 - (2\omega)^2\mathbf m_1\right]\theta_2 = \mathbf f_3 \\ &\quad \vdots \end{align}
Then I can turn the above into a matrix:
\begin{equation} \begin{bmatrix} \mathbf A_{1,1} \\ \left[\mathbf g_1 - (\omega)^2\mathbf m_1\right] & \mathbf A_{2,2} \\ \left[\mathbf g_2 - (\omega)^2\mathbf m_2\right]& \left[\mathbf g_1 - (2\omega)^2\mathbf m_1\right] & \mathbf A_{3,3} \\ \vdots & \vdots & \vdots & \ddots \\ \end{bmatrix} % \begin{bmatrix} \theta_1 \\ \theta_2 \\ \theta_3 \\ \vdots \end{bmatrix}= % \begin{bmatrix} \mathbf f_1 \\ \mathbf f_2 \\ \mathbf f_3 \\ \vdots \end{bmatrix} \tag{4} \end{equation}
where for rows with index $r$ and columns with index $c$ we have:
\begin{gather} \mathbf A_{r,r=c}=ir\omega\mathbf C - (r\omega)^2\mathbf M_v + \mathbf K \\ A_{r, c<r} = \mathbf g_{r-c} - (c\omega)^2\mathbf m_{r-c} \end{gather}
where the above is a diagonal matrix since functions $\mathbf M_v$ and $\mathbf G$ are also diagonal.
Having solved Eq. 4 to obtain the coefficients $\theta_n$ for the positive part of the sum, the conversion back to the time domain can be done as follows:
$$\mathbf\Theta(t) = \sum^{N}_{n=1} 2|\theta_n|\cos(n\omega t - \phi_n)$$
where $\phi_n = \text{arg}(\bar{\theta}_n)$, where $\bar{\theta}_n$ is the conjugate of $\theta_n$.