Task 5: A more sophisticated IMEX solver for the Cahn-Hilliard equation

Before we finally run some more interesting simulation of the actual Cahn-Hilliard equation with realistic initial conditions (and no source term!), we want to improve the previous IMEX solver for the Cahn-Hilliard equation and implement a more sophisticated solver for the Cahn-Hilliard equation which can be shown to be second-order accurate in time.

Again, the scheme is best explained by looking at a general ODE system of the form

\[ U_t = \mathbf{L}U + \mathbf{N}(U), \quad U(0) = U_0, \]

Now the Song scheme from [Song, 2016] is a 3-stage IMEX Runge-Kuttta scheme given by

(89)\[\begin{align} U^{(1)} &= U^n + \tau\bigl( \mathbf{L}U^{(1)} + \mathbf{N}(U^n) \bigr), \\ U^{(2)} &= \alpha_{10} U^n + \alpha_{11} U^{(1)} + \beta_1 \tau\bigl( \mathbf{L}U^{(2)} + \mathbf{N}(U^{(1)}) \bigr), \\ U^{n+1} &= \alpha_{20} U^n + \alpha_{21} U^{(1)} + \alpha_{22} U^{(2)} +\beta_2 \tau\bigl( \mathbf{L}U^{n+1} + \mathbf{N}(U^{(2)}) \bigr), \end{align}\]

Theoretical tasks

1) ~~Show that this~~ This scheme has consistency order \(p=2\) if the coefficients \(\beta_i, \alpha_{ij}\) satisfies the following order conditions:

\[\begin{split} \left\{ \begin{aligned} \alpha_{10} + \alpha_{11} &= 1, \\ \alpha_{20} + \alpha_{21} + \alpha_{22} &= 1, \\ \alpha_{21} + \alpha_{22} \alpha_{11} + \alpha_{22} \beta_1 + \beta_2 &= 1, \\ \alpha_{21} + \alpha_{22} \alpha_{11} + \alpha_{22} \alpha_{11} \beta_1 + \alpha_{22} \beta_1^2 + \alpha_{21} \beta_2 + \alpha_{22} \alpha_{11} \beta_2 + \alpha_{22} \beta_1 \beta_2 + \beta_2^2 &= \frac{1}{2}, \\ \alpha_{22} \beta_1 + \alpha_{11} \beta_2 + \beta_1 \beta_2 &= \frac{1}{2}. \end{aligned} \right. \end{split}\]

~~To arrive at these equations, employ a similar Taylor-expansion technique as shown in one of the homework exercises and in the derivation of the consistency order of Heun’s method and the improved Euler methods.~~

This set of equations cannot be solved uniquely. Here, we consider the following coefficients which solve the above equations:

(90)\[\begin{align} \alpha_{10} &= \frac{3}{2}, \quad &\alpha_{11} &= -\frac{1}{2}, \quad &\alpha_{20} &= 0, \quad &\alpha_{21} &= 0, \quad &\alpha_{22} &= 1, \quad &\beta_1 &= \frac{1}{2}, \quad &\beta_2 &= 1. \\ \alpha_{10} &= 2, \quad &\alpha_{11} &= -1, \quad &\alpha_{20} &= \frac{1}{2}, \quad &\alpha_{21} &= 0, \quad &\alpha_{22} &= \frac{1}{2}, \quad &\beta_1 &= 1, \quad &\beta_2 &= 1. \\ \alpha_{10} &= 2, \quad &\alpha_{11} &= -1, \quad &\alpha_{20} &= 0, \quad &\alpha_{21} &= \frac{1}{2}, \quad &\alpha_{22} &= \frac{1}{2}, \quad &\beta_1 &= 1, \quad &\beta_2 &= \frac{1}{2}. \\ \alpha_{10} &= \frac{5}{2}, \quad &\alpha_{11} &= -\frac{3}{2}, \quad &\alpha_{20} &= \frac{2}{3}, \quad &\alpha_{21} &= 0, \quad &\alpha_{22} &= \frac{1}{3}, \quad &\beta_1 &= \frac{3}{2}, \quad &\beta_2 &= 1. \end{align}\]

2) Apply this scheme now to the the Cahn-Hilliard equation formulated in the Fourier space. Here, in each of the above stages, you should apply the same convex splitting as in the previous task. As always, before you start implementing the scheme, provide a brief mathematical description of the resulting numerical scheme, in particular, describe how a new solution is computed in the Fourier space from the previous solution for each time step.

Computational tasks

3) Rerun the EOC study from the previous task for the new IMEX scheme for all 4 sets of coefficients, compare the results against each other and with the previous IMEX scheme. Which of the 4 Song coefficients sets will you favor for the remaining project? Here, you can focus on the case \(\kappa = 0.01\).

Important note: To take into account an inhomogeneous right-hand side \(\mathbf{G}\), you can implement the following slightly modified scheme:

(91)\[\begin{align} U^{(1)} &= U^n + \tau\bigl( \mathbf{L}U^{(1)} + \mathbf{N}(U^n) + \mathbf{G}^{n+1/2} \bigr), \\ U^{(2)} &= \alpha_{10} U^n + \alpha_{11} U^{(1)} + \beta_1 \tau\bigl( \mathbf{L}U^{(2)} + \mathbf{N}(U^{(1)}) + \mathbf{G}^{n+1/2} \bigr), \\ U^{n+1} &= \alpha_{20} U^n + \alpha_{21} U^{(1)} + \alpha_{22} U^{(2)} +\beta_2 \tau\bigl( \mathbf{L}U^{n+1} + \mathbf{N}(U^{(2)}) + \mathbf{G}^{n+1/2} \bigr), \end{align}\]

where \(\mathbf{G}^{n+1/2} = \mathbf{G}(t_n + \tau/2)\).