Task 3: A spectral solver for the transient biharmonic equation#
As the third step towards a solver for the Cahn-Hilliard equation, you are now ask to implement, test and verify a solver for the time-dependent biharmonic equation
As before, 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, i.e. the same considerations as for the biharmonic equation apply. Moreover, the problem is supplemented with initial conditions \(u(0,x,y) = u_0(x,y)\), which we also assume to satisfy the periodic boundary conditions.
This problems resembles very closely the heat equation, the only difference being that the spatial differential operator is now given by the biharmonic operator, and not \(-\Delta\) as in the heat equation. Consequently, we take a similar approach as we took for the heat equation in the lecture, and combine a Fourier spectral method in space with a one-step time-stepping methods in time. The one-step method we ask you to use here is called the \(\theta\)-method.
The \(\theta\)-method#
Let’s forget for the moment that we want to solve the biharmonic equation, and consider a general linear ODE of the form
The \(\theta\)-method is defined as follows. As usual, we first discretize the time interval \([0,T]\) into \(N\) equidistant time steps of size \(\tau = T/N\). The algorithm is started from the initial condition \(U^0\). at time \(t_0\). Then, given the solution \(U^n\) at time \(t_n\), the new solution \(U^{n+1}\) at time step \(t_{n+1} = t_n + \tau\) is computed as
Theoretical tasks#
1) Rewrite the \(\theta\)-method as a 2-stage Runge-Kutta method, derive the corresponding Butcher table, and discuss the consistency order of the \(\theta\)-method for various values of \(\theta\). For which values does the \(\theta\)-method reduces to other well-known time-stepping methods?
Hint: Have a look at the order conditions for Runge-Kutta methods in the lecture notes.
2) Next, determine the stability function \(r_{\theta}(z)\) of the \(\theta\)-method and plot the stability region of the \(\theta\)-method in the complex plane for \(\theta = 0, 0.25, 0.498, 0.5, 0.502, 0.75, 1\). For which values of \(\theta\) does the \(\theta\)-method seem to be A-stable? What do you conjecture for a general \(\theta\)? Do you have an idea of how the stability region of the \(\theta\)-method looks in general/depends on \(\theta\)?
3) Finally, we ask you to put your conjecture on solid mathematical grounds. More precisely, determine mathematically the stability region of the \(\theta\)-method and how it depends on the value of \(\theta\).
Hint: The border between the stable and unstable region is given by \(\partial S_{\theta} = \{z \in \mathbb{C} : |r_{\theta}(z)| = 1 \}\) and this can be transformed into a simple equation for a circle in the complex plane. How does the center and radius of this circle depend on \(\theta\)?
Computational tasks#
After this theoretical warm-up, we now turn to the implementation of the spectral solver for the transient biharmonic equation. Please implement a solver for the transient biharmonic equation combining the Fourier spectral method in space with the \(\theta\)-method in time.
1) Before you start implementing, please provide a brief mathematical description of the resulting numerical scheme, describing how a new solution is computed from the previous solution for each time step.
For the implementation, the following specifications for the solver interface be met:
The solvers should be implemented as a generator function using the yield statement to return the
discrete Fourier transform of the solution at each time step, and the current time.
The generator function should have the following signature:
def transient_biharmonic_solver(*, kappa,
X, Y, U0,
t0, T, Nt,
theta=0.5,
g=None):
"""
Solve the transient biharmonic equation using the theta method.
Parameters:
-----------
kappa (float): Diffusion coefficient.
X (ndarray): 2D array of x-coordinates.
Y (ndarray): 2D array of y-coordinates.
U0 (ndarray): Initial condition array.
t0 (float): Initial time.
T (float): Final time.
Nt (int): Number of time steps.
g (callable, optional): Source term function g(X, Y, t). Defaults to None.
Yields:
-------
tuple: A tuple containing the discrete Fourier transform of U at t, and the current time t.
"""
# Prepare relevant data for Fourier transform
...
...
# Compute DFT of the initial value
...
# Add your time-stepping loop here
# Time stepping
t = t0
dt = (T-t0)/Nt
# For convenience when plotting, computing errors, etc.,
# return the initial solution and initial time.
yield U_hat, t
while t < T-dt/2:
# Solve for next time step and update time
...
...
yield Uhat, t
which then can be used to solve the transient biharmonic as follows:
# Set up the problem
...
...
solver = transient_biharmonic_solver(kappa=kappa,
X=X, Y=Y, U0=U0,
t0=t0, T=T, Nt=Nt,
theta=theta,
g=None)
for Uhat, t in solver:
# Do something with the solution
...
...
To verify your implementation and assess the stability and accuracy of the solver, we ask you perform the following tasks:
2) First, study the convergence order of the time-stepping scheme for \(\theta \in \{0, 0.5, 1\}\). Use the function \(u_{\mathrm{ex}}(x,y,t) = \sin(x)\cos(y)\exp(-\lambda \kappa t)\) to manufacture a solution for the transient biharmonic equation. Choose \(\lambda\) such that the manufactured solution leads to a homogeneous source term \(g(x,y,t) = 0\).
Set \(\Omega = [-\pi,\pi)^2\) and \(\kappa = 1\) and choose \(N = N_x = N_y = 20\) sampling points/subintervals in each space direction. Furthermore, set \(t_0 = 0\), \(T = 1\). Now solve the problem successively for \(N_t = 10, 20, 40, 80, 160, 320, 640\) time steps with equidistant time steps \(\tau = T/N_t\). For each run calculate the error in the so-called \(L^{\infty}L^{\infty}\) norm defined by
and compute the EOC with respect to the time step size \(\tau\) (or equivalently the number of time steps \(N_t\)).
Display your results in a table showing the number of steps, resulting error and experimentally observed convergence order for each refinement in time. Do this for \(\theta = 0, 0.5, 1\), starting with \(\theta = 1\), then \(\theta = 0.5\) and finally \(\theta = 0\).
Discuss the results and relate them to the theoretical results you derived in the first part of the task. Comment on the suitability of the \(\theta=0\)-method for the transient biharmonic equation. In particular, explain why for \(\theta = 0\) the scheme fails by deriving the resulting CLF condition for this case.
Finally, based on the derived CFL condition find the minimal number \(N_{CLF}\) of time steps \(N_t\) for which the scheme should be stable and repeat the convergence study for \(\theta = 0\) with \(N_t = 0.5 N_{CLF}, N_{CFL}, 2 N_{CFL}, 4 N_{CFL}\).
3) Finally, to prepare yourself for computing manufactured solutions for the Cahn-Hilliard equation in the next task, we ask you to rerun the convergence study for the transient biharmonic equation for \(\theta = 0.5, 1\) with the manufactured solution
resulting in a highly non-trivial source term \(g(x,y,t)\). Use \(N_t = 10, 20, 40, 80, 160, 320, 640\) as before for the EOC study.
Hint: Make your live easy and use the sympy to manufacture solutions by computing the resulting source terms. See the example from the lecture notes.