Task 2: A spectral solver for the biharmonic equation

Next, following the example of the Poisson solver from the lecture, we ask you to implement a spectral solver for the biharmonic equation supplemented with a \(0\)-order term in \(2d\).

The resulting equation reads

\[ \Delta^2 u + c u = f \quad \text{in } \Omega, \]

on a rectangular domain \(\Omega = [0,L_x)\times [0,L_y)\), where \(\Delta^2\) is the biharmonic operator and \(c \geqslant 0\) is a non-negative constant.

The problem is to be solved on a given rectangular domain \(\Omega = [0,L_x) \times [0,L_y) \subset \mathbb{R}^2\) with periodic boundary conditions. More precisely, we assume that all derivatives up to order 3 are periodic. Different side lengths \(L_x\) and \(L_y\) are allowed and should be incorporated in the implementation.

Now implement a spectral solver for the biharmonic equation using the fast Fourier transform. The solver should be implemented in a function of the form

def biharmonic_solver(X, Y, F, c, mean=0.0):
    """
    Solve the biharmonic equation in 2D using the spectral method.

    Parameters:
        X (ndarray): 2D array of x-coordinates.
        Y (ndarray): 2D array of y-coordinates.
        F (ndarray): 2D array representing the right-hand side of the biharmonic equation.
        c (float): Constant coefficient in the biharmonic equation.
        mean (float, optional): Desired mean value of the solution in case c = 0. Default is 0.0.

    Returns:
        U (ndarray): 2D array representing the solution to the biharmonic equation.
    """
    pass # Add your code here

Your solver should also be able to gracefully handle the case \(c=0\) by prescribing a (user-defined) mean value for the solution. This will come in handy when you are asked to test your solver using manufactured solutions which do not necessarily have a 0 mean value.

To verify the correctness of your implementation, run a number convergence studies using the manufactured solutions from the following 2 test cases, each of them posed on the rectangular domain \(\Omega = [0,2\pi) \times [0,4\pi)\):

• \(u(x,y) = \sin(8(x-1))\cos(4y)\), \(c=1\), and \(N_x = 4, 8, 15, 16, 20, 32\) and \(N_y = 2N_x\).

• \(u(x,y) = \exp(\sin^2(x) + \cos(2y))\), \(c=0\), \(N_x = 4 + 4k\), \(k=0,1\ldots,9\) and \(N_y = 2N_x\).

Remember to compute the corresponding right-hand side \(f\) for the manufactured solutions. For each series, compute the experimental order of convergence (EOC) with respect to the maximum norm over the grid points report them in a table.

Provide also a surface plot of both the exact solution, the numerical solution and the error function for \(N_x = 15, 16\) in the first series and \(N_x = 32\) in the second series.

Discuss the results.

Hints:

• To compute the right-hand side \(f\) for the manufactured solutions, it might be helpful to start using the sympy module, see lecture notes and the tutorial. This will definitely pay-off later when you are asked to manufacture solution for the Cahn-Hilliard equation!

• Recall that numerically the check \(c=0\) is not trivial. As a rule of thumb, you can check e.g. whether \(c\) is smaller the smallest non-zero quartic frequency in the Fourier space, i.e., \(c < \mathrm{tol} \cdot k_{\text{min}}^4\), where \(k_{\text{min}}\) is the smallest non-zero frequency in the Fourier space.