Solving Cahn-Hilliard equation using semi-implicit Fourier spectral methods

Question

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

Answer 1

Based on the references you provided, I assume your initial conditions are random numbers (see second reference, Sect. IV)? If so, it's probably no surprise that you do not see higher-order convergence in time, at least not from the beginning. Your initial conditions are not smooth and it may take some time to get to a smooth solution/state. There, however, higher-order convergence should be seen. I recommend starting with a smooth and periodic initial condition to check convergence order.

In case your initial conditions are different (i.e. smooth and periodic), please specify them and share the code, ideally as a minimal working example.