Cannot specify full boundary conditions in biharmonic equations without periodic boundary conditions

[tranqui]

I am trying to solve 4th order ODEs containing a biharmonic operator. The documentation describes solving the Kuramoto-Sivashinsky equation and the Cahn-Hilliard equation.

However, there is no option to specify enough boundary conditions to solve a 4th order ODE. For a linear ODE we need 4 and this is generally true for nonlinear ODEs also. With Neumann BCs I expect to be able to specify values and derivatives at the same point, with e.g. something like for a 2d system:
bcs = [{'value': -1, 'derivative': 0}, {'value': 1, 'derivative': 0}]
However, this always throws an error saying this feature is not supported. So I am curious how 4th order ODEs are listed as examples. Ultimately I'm interested in something similar to the Cahn-Hilliard equation so look to that as an example.

Looking at the implementation of the Cahn-Hilliard equation, it looks like this would only work with periodic boundary conditions because this implicitly sets values and all derivatives. That being said, I am never able to get convergent results for e.g. the binodal so I suspect either I am using py-pde wrong or there is a problem with the underlying solver. Below is a minimal working example that illustrates this where I try to obtain the binodal. I start with two interfaces around x=+-0.5 to respect periodic boundaries at x=+-1. The solver quickly fails with a "Field was not finite" error:

from pde import CahnHilliardPDE, ScalarField, CartesianGrid
grid = CartesianGrid([[-1, 1], [-1, 1]], [1000,50])
state = ScalarField.from_expression(grid, 'tanh(cos(pi*x)/0.1)')
eq = CahnHilliardPDE()
result = eq.solve(state, 1, dt=1e-6, method='implicit')
result.plot()

My expected behaviour is that it should converge to tighten up the interfaces with the correct interface width (rather than the dummy value of 0.1 I gave as an IC).

Is this a bug, or am I using it wrong? If the latter could you please advise me on how I should use the Cahn Hilliard (and similar) solvers? Thank you.

[david-zwicker]

You're asking multiple questions instead of describing a particular issue with the code, so I converted it into a discussion item. I'll try to address the questions I could find:

• A biharmonic operator is basically two nested Laplace operators. Each Laplace operator requires a boundary condition at each boundary point in py-pde. The CahnHilliardPDE class allows setting these conditions using the bc_c and bc_mu flag. When using the PDE class, you can set specific boundary conditions for specific operators using the bc_ops argument; see the documentation.
• I don't know what you tried to achieve when setting the boundary condition bcs = [{'value': -1, 'derivative': 0}, {'value': 1, 'derivative': 0}], but I'm not surprised that you got an error.
• Your code for the Cahn-Hilliard equation will likely not converge because you discretization length is orders of magnitude smaller than the interface length. A working example reads:

from pde import CahnHilliardPDE, ScalarField, CartesianGrid
grid = CartesianGrid([[-10, 10], [-10, 10]], 30)
state = ScalarField.from_expression(grid, 'tanh(cos(0.1 * pi*x)/0.1)')
eq = CahnHilliardPDE()
result = eq.solve(state, 10, dt=1e-6, adaptive=True)
result.plot()

I hope this gets you on the right track!

[tranqui]

Many thanks, that's a very helpful document! I think I see the problem. Normally in order to go to fourth order one has to introduce 2 ghost points on each boundary, so indices $\{0, 1, N-2, N-1\}$ would become the ghost points. Then when you evaluate the laplacian twice (or otherwise use a next-nearest neighbour stencil) you can still apply two BCs on the edge points. So I guess that would require modifying the grid class (if that's where ghost points are encoded?).

If the architecture is such that adding additional ghost points would be unwieldy, perhaps a hacky solution would be to switch to a non-central stencil at the boundaries when additional BCs are to be applied. Either a forward/backward or a mixed stencil with one ghost point would do the trick. That would still require next-nearest neighbours (in the inward direction) though to achieve second-order accuracy.

Equation (11) of the document linked above puzzles me. It looks like you're only applying first-order accurate first derivatives at the boundary, and not strictly taking the boundary value $f(0) = y_1$ but rather the mid-point between it and the virtual point. Why not simply continue the second-order stencil?

[tranqui]

Right, you apply BCs at the operator level so you have to precalculate the ghost points. So they're only ever known to first-order. Introducing additional ghost points would then potentially improve the overall order of the solver then?

[tranqui]

Though you could use $f(0) = \frac{2y_0 + y_{-1} + y_1}{4} + \mathcal{O}(\Delta x^2)$ and $f'(0) = \frac{y_1 - y_{-1}}{2\Delta x} + \mathcal{O}(\Delta x^2)$. It's straightforward to find similar expressions for higher-order ghost points - I can provide you with expressions/code if that would be helpful?

[david-zwicker]

Note that we do not have support points on the boundary but rather the boundary is half way between support points; see the documentation. This has many advantages and allows to express values and derivatives using Eq. 11 in the PDF document. I agree that a fourth-order derivative could in principle be implemented with a five-point stencil, but since this requires two virtual points in each direction, we didn't do this. Instead, we implement the biharmonic operator by applying the Laplacian twice, which lets us easily impose boundary conditions for u and u''. There might be a way to calculate the virtual point for u'' from the boundary condition on u', but this would probably also involve knowledge about the actual differential equation.

[tranqui]

I think I understand a bit better now. I doubt there is a clean way to implement this kind of boundary condition purely at the level of the u'', as it affects the other differential operators equally. As you say, you'd need some knowledge of the underlying PDE to "undo" ignorance of this BC in the other operators. I can't see a way of doing it without switching to a higher-order stencil unfortunately, which may be too unwieldy to be worthwhile.