Consider the integral
$$2\pi T = -\frac{1}{2} \int_0^{2\pi} A\frac{\partial B}{\partial \xi} \ d\xi$$
where $$A= \sum_{-N}^N i \ sign(n) \ B_n e^{-in\xi}; \quad \quad B= \sum_{-N}^N \ B_n e^{-in\xi}.$$ These expansions imply
$$2\pi T = \sum |k|B_kB_{-k}$$.
Furthermore, $$B_k = \frac{i}{k}\left(-\dot{Y}_k + \sum_{l+j=k}\dot{Y}_lY_j\ j (sign(l)-sign(j))\right),$$
with dots indicating temporal derivatives.
Now, the canonical variables of the system end up being the complex valued variables $Y_{\pm 1},Y_{\pm 2},..,Y_{\pm N}$, and so the Euler-Lagrange equations involve terms like
$$ (1) \quad \quad \frac{\partial T}{\partial Y_m}; \quad \quad \frac{d}{dt}\left(\frac{\partial T}{\partial \dot{Y}_m}\right).$$
Therefore, to derive the governing equations, I must explicitly expand terms like these and write them as some function $F(Y_{\pm 1},Y_{\pm 2},..., \dot{Y}_{\pm 1},...)$
Now, the way I'm currently computing this term is to note
$$T = \sum_{k,l} Q_{kl} \dot{Y}_k\dot{Y}_l,$$
where $$Q_{kl} = \frac{1}{4}\sum_{m=1}^N P_{mk}P_{-ml}$$
and $$P_{mn} = (sign(m-n)-sign(n))Y_{m-n} -\frac{1}{2} Y_mY_{-n}.$$
To set up my governing equations via the Euler-Lagrange equations, I need to compute terms like the ones mentioned above in (1). They both lead to the need to calculate
$$\frac{\partial Q_{kl}}{\partial Y_n} = \sum_{m=1}^N \left(\frac{\partial P_{mk}}{\partial Y_n}P_{-ml}+P_{mk}\frac{\partial P_{-ml}}{\partial Y_n}\right).$$
I need to compute this for all n, so I end up needing to perform a series of 2N+1, $((2N+1),N)\times (N,2N+1))$ matrix multiplications. For large $N$, this computation is absolutely slaughtering my run times, and generally take up +80% of the run time for each time step (depending on the size of N). This is my problem.
The way I've been trying to overcome this is to focus my attention on using the best optimized BLAS libraries to compute the matrix multiplications, and exploiting symmetries of the matrices to try to reduce computation time. However, I think I've milked this line of reasoning for all it's worth without getting too deep into the weeds, so I would like to atleast think about alternatives.
I am very new to numerical methods, and might be completely off base, but it seems like there should be a way to manipulate the FFT here, to avoid doing these matrix multiplications as these Cauchy products are similar to convolutions (ie the pseudo spectral method).
However, I don't see how one can explicitly do this, especially as I need to differential these products with respect to the dependent variables of the system.
If anyone has any insight, it's greatly appreciated.
