So, I have written both a C and python code to solve the 2D Cahn-Hilliard equation: \begin{equation} \frac{\partial c}{\partial t} = \nabla^2\left(c^3 - c - \kappa\nabla^2c\right) \end{equation} Taking its Fourier transform yields the following ODE: \begin{equation} \frac{\partial \hat c}{\partial t} = \underbrace{-k^2\{c^3-c\}_\mathbf{k}}_{F\left(\hat c\right)} \overbrace{- \kappa k^4\hat c}^{G\left(\hat c\right)} \end{equation} Where G is linear and "fast-varying" (high eigenvalues) and F "slow-varying" (low eigenvalues) and $\mathbf k = \sqrt{k_x^2 + k_y^2}$ the Fourier "wavelength". To overcome the highly-constraining 4th-order term in terms of timestep size, I decided to use the following second order implicit-explicit "IMEX" scheme where $\hat f \equiv \{c^3 - c\}_\mathbf k$: \begin{equation} \left(3+2\Delta t\kappa k^4\right)\hat c^{n+1} = 4\hat c^n - \hat c^{n-1} - 2\Delta t k^2\left(2\hat f^n - \hat f^{n-1}\right) \end{equation} Where the iterations are started using the first order IMEX-BDF scheme. Now, this works perfectly well and exhibits second-order error convergence as expected. I then wanted to try and use the fourth-order IMEX-BDF scheme to obtain even better accuracy at basically no extra CPU cost (it requires more memory to store the previous steps, but memory is not the issue in this case): \begin{equation} \left(25+12\Delta t\kappa k^4\right)\hat c^{n+1} = 48\hat c^n - 36\hat c^{n-1} + 16\hat c^{n-2} - 3\hat c^{n-3} - 12\Delta t k^2\left(4\hat f^n - 6\hat f^{n-1} + 4\hat f^{n-2} - \hat f^{n-3}\right) \end{equation} However, surprisingly to me, this scheme yields a slightly worst second order error convergence behavior than the second order scheme (slightly worst in that the error constant looks to be bigger for the IMEX-BDF4 scheme). I have observed the same behavior for the IMEX-BDF3 scheme (second order error convergence instead of third; slightly worst error constant than IMEX-BDF2 - but a bit better than IMEX-BDF4).
To validate the convergence order, I took as reference solution a fourth-order classical Runge-Kutta scheme with very low $\Delta t$.
Thus my question is the following: is there any reason why the 4th- and 3rd-order IMEX-BDF schemes would be somehow limited to second order in the particular use case of the Cahn-Hilliard equation ? Or is it me who implemented something the wrong way ? Since I'm using Fourier spectral methods for physical differentiation, I'm assuming the physical discretisation is not a limiting factor in terms of simulation precision since I'm using a fairly high number of mesh points - 256x256.
My python code w/ IMEX2, IMEX4 & RK4: pastebin
References:
Willem Hundsdorfera, Steven J. Ruuth, IMEX extensions of linear multistep methods with general monotonicity and boundedness properties
Jingzhi Zhu, Long-Qing Chen, Jie Shen, and Veena Tikare, Coarsening kinetics from a variable-mobility Cahn-Hilliard equation: Application of a semi-implicit Fourier spectral method