The Problem
I currently have the following system of ODEs which I have managed to solve successfully using the ode23tb solver in MATLAB.
$$(\mathbf{J_c}+\mathbf{J_v}(t))\mathbf{\ddot{\Theta}} + \mathbf{C\dot{\Theta}}+\mathbf{K\Theta} = \mathbf{F}(t) \tag{1}$$
However, using ode23tb (or similar solvers) takes time, and the solution has to go through a transient stage before it settles into a repeatable periodic pattern. Hence I have attempted to use Fourier series which might bypass the above problems.
My Attempt
- I attempt to solve the problem iteratively where for the first iteration with superscript $^{(1)}$ I take the Fourier transform of the RHS of Eq.(1), i.e. $\mathbf{F}^{(1)}$:
$$\mathbf{F}(t) = \mathbf{F}^{(1)} = \sum^{+\infty}_{n=-\infty} \mathbf{f}^{(1)}_n e^{in\omega t}$$
- I remove the time dependence of the $\mathbf{J}_v$ term by taking its average $\mathbf{\bar{J}}_v$, and perform the Fourier transform on Eq.(1) to get a linear set of equations:
$$(-n^2 \omega^2 (\mathbf{J}_c+\mathbf{\bar{J}}_v) + i n \omega \mathbf{C} + \mathbf{K})^{-1} \mathbf{f}^{(1)}_n = \mathbf{c}^{(1)}_n \tag{2}$$
- Having solved Eq.(2), I can now use the fact that $ \ddot{\mathbf{\Theta}}^{(1)} = \sum -(n\omega)^2 \mathbf{c}^{(1)}_n e^{in\omega t}$ to re-calculate the excitation force on the RHS of Eq.(1). So for my second iteration Eq.(1) becomes: $$\mathbf{J_c}\mathbf{\ddot{\Theta}} + \mathbf{C\dot{\Theta}}+\mathbf{K\Theta} = \mathbf{F}(t)-\mathbf{J_v}(t)\mathbf{\ddot{\Theta}}^{(1)} \tag{2}$$ Therefore I repeat steps 1 and 2 by setting as $\mathbf{F}^{(2)}$ the entire RHS of Eq.(2) above.
My hope is that after a couple iterations, the values of subsequent solutions of Eq.(2) will converge however the solution amplitudes become higher and higher.
I have also compared the acceleration with the one obtained from ode23tb for the first iteration, and they are nearly identical, which shows that at least the first iteration is performed correctly. I sense its my method that is at fault here.