I'm trying to model some two-dimensional waves, and am unsure how to combine my boundary conditions with spectral methods. The PDE I'd like to explore resembles the equations for damped, driven bending-mode waves in a thin disk or beam:

$$\ddot w(\boldsymbol x,t) = -c \nabla^2 \nabla^2 w(\boldsymbol x,t) - d \dot w(\boldsymbol x,t) + u(\boldsymbol x, t),$$

where $w(\boldsymbol x,t)$ is the spatiotemporal response of the system (analogous to vertical displacement of a horizontally oriented thin plate), $\boldsymbol x = \{x_1,x_2\}$ are 2D spatial coordiantes, $u(\boldsymbol x, t) = u(\boldsymbol x) u(t)$ is a spatiotemporally separable external input, and $- d \dot w(\boldsymbol x,t)$ is a damping/dissipation term with unknown $d$.

I've settled on this equation for my model based on an empirical dispersion relationship seen in experimental data; There may be other equivalent models. For example, the Schrödinger equation has a similar dispersion relationship and I'd be equally happy with methods employing that if its any easier.

The boundary I'd like to use looks like this:

boundary conditions

The actual region I want to simulate isn't circular (or defined by any simple function), this was just easier to draw.

This is what I've done so far:

  • I define a square 2D grid with resolution sufficiently smaller than the wavelengths of interest.
  • I mask out a sub-domain of this grid with the irregular shape I'd like to model.
  • I construct a graph adjacency matrix for the pixels/cells in this grid.
  • I calculate the eigenvectors and eigenvalues of this matrix, and keep the top 16 components as a low-rank basis on which to approximate my solutions.
  • These 16 components represent the top 16 eigenmodes of the Laplacian operator $\nabla^2$ with Dirichlet boundary conditions set to zero.

From here, numerical solutions are relatively straightforward. I prefer to think of the functions $w$ and $u$ as vectors in a Hilbert space, and write the spatial and time dependances as subscript indecies:

$$\ddot {\boldsymbol w}_{\boldsymbol x,t} = - a {\nabla^2 \nabla^2} {{\boldsymbol w}_{\boldsymbol x,t}} -d \dot {\boldsymbol w}_{\boldsymbol x,t} + {\boldsymbol u}_{\boldsymbol x} {\boldsymbol u}_t. $$

Take the Fourier transform in time , and change to the eigenbasis of the Laplacian operator in space. Denote temporal angular frequency as $\nu$. Denote spatial eigenmode number as $k$, with associated eigenvector $\boldsymbol v_{k}$ and eigenvalue $-\lambda_k$. We'll denote the thusly transformed components as $\boldsymbol\xi_{k,\nu}$. Note that the Fourier transform of the first and second (time) derivatives are $i\nu$ and $-\nu^2$, respectively. In this basis, the PDE becomes:

$$-\nu^2 {\boldsymbol\xi}_{k,\nu} = -a \lambda_k^2 \boldsymbol\xi_{k,\nu} - i\nu d {\boldsymbol\xi}_{k,\nu} + \hat{\boldsymbol u}_k \tilde{\boldsymbol u}_\nu,$$

where $\hat{\boldsymbol u}_k$ is the projection of $\boldsymbol u_{\boldsymbol x}$ onto the $k^{\text{th}}$ spatial mode and $\tilde{\boldsymbol u}_\nu$ is the (temporal) Fourier transform of ${\boldsymbol u}_t$. Collect terms and solve for ${\boldsymbol\xi}_{k,\nu}$:

$${\boldsymbol\xi}_{k,\nu} = \frac {\hat{\boldsymbol u}_k \tilde{\boldsymbol u}_\nu} {a \lambda_k^2 + i\nu d -\nu^2}$$

The spatiotemporal solution can then be obtained by applying the inverse Fourier transform to ${\boldsymbol\xi}_{k,\nu}$ and projecting back onto our 16 spatial components. This approach is much faster than simulating and waiting for said simulation to reach steady-state, and lets me explore parameter variations efficiently.

The issue, however is: how can one apply the same approach using the mixed ("Robin"?) boundary conditions described earlier? I suspect these boundary conditions should change the eigendecomposition, but I'm not sure how to construct/describe an operator encoding them?

Many thanks in advance! Disclaimer: I don't normally do this sort of thing, so apologies if I've mis-used some of the terminology.


1 Answer 1


For this particular finite-element approach, Dirichlet (zero) boundary conditions can be implemented by zeroing-out the Laplacian operator at boundaries; Neumann (zero derivative) boundary conditions can be implemented using the graph Laplacian of the simulated sub-domain (implemented such that out-of-bounds points aren't counted as neighbors).

To clarify, consider an example in 1D:

  • In the interior, the discrete Laplacian is $\{1,-2,1\}$
  • At an e.g. right edge for Dirichlet boundary conditions, the discrete Laplacian is $\{1,-2,0\}$
  • For Neumann boundary conditions, the discrete Laplacian at this boundary is $\{1,-1,0\}$.

The boundary conditions for a 2D finite-element approach are analogous.

  • You can interpolate between reflective (Dirichlet) and absorbing (Neumann) boundary conditions by taking a linear combination of their respective discrete Laplacian operators at boundaries.
  • You can implement a gradual transition between these two boundary conditions by smoothly changing this interpolation weights along the edge of your simulation.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.