# Applying spectral method to a damped, driven 2D bending-mode wave equation on an irregular domain with heterogeneous boundary conditions

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: 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.

• 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\}$$.