The following is the Discrete Spectral Vorticity Evolution PDE for incompressible flow: $$ \frac{\partial \Omega_{pq}}{\partial t} = \nu \left( \frac{\partial^2 \Omega_{pq}}{\partial x^2} + \frac{\partial^2 \Omega_{pq}}{\partial y^2} \right) - \mathrm{fft2} \left( u_{ij} \frac{\partial \omega_{ij}}{\partial x} + v_{ij} \frac{\partial \omega_{ij}}{\partial y} \right) $$ After significant manipulation, we can obtain the following equivalent form, where $\Xi_{pq} = \exp \left[ \nu \left( k_p^2 + k_q^2 \right) t \right]$ is an integrating factor $$ \frac{\partial}{\partial t} \Big[ \Xi_{pq} \Omega_{pq} \Big] = -\underline{i} \Xi_{pq} \Big[ k_p \cdot \mathrm{fft2} \left(u_{ij} \omega_{ij} \right) + k_q \cdot \mathrm{fft2} \left( v_{ij} \omega_{ij} \right) \Big] $$ This secondary form is advantageous for numerical integration because we no longer need to consider both first (time derivative) and second-order (viscous term) stability, but brings with it a new problem.
The wavenumbers $k_p$ and $k_q$ are defined as follows (assuming a domain from $[0, 2 \pi]$, with $M$ points not including at $2 \pi$): $$ k_p = \begin{cases} p - 1, \quad \quad \quad \, p \leq \mathrm{floor} \left( \frac{M}{2} \right) + 1 \\ -M + p - 1, \quad p > \mathrm{floor} \left( \frac{M}{2} \right) + 1 \end{cases}$$ Because the wavenumbers are squared, our integrating factor $\Xi_{pq}$ is constantly increasing with time, and even though our viscosity constant $\nu$'s small size prevents an immediate blow-up, eventually we do start having overflow problems as time increases.
Is there any way to prevent this while still only maintaining first-order stability?
One idea I had was to use some type of logarithm-exponent trick, but that ran out of steam when I didn't know what to do with the LHS.
