# Solving a biharmonic eigenvalue Problem

I am interested in solving the following biharmonic eigenvalue problem.

$$\begin{array}{cccc} & \Delta ^2 \Psi (x,y) = \lambda \Psi (x,y), & - \frac{\pi}{2} \le x \le \frac{\pi}{2} & -\frac{\pi}{2} \le y \le \frac{\pi}{2} \\ & x = \phantom{-}\frac{\pi}{2} & \Psi = 0 & \dfrac{\partial \Psi }{\partial x} = 0 \\ & x = - \frac{\pi}{2} & \Psi = 0 & \dfrac{\partial \Psi }{\partial x} = 0 \\ & y = \phantom{-}\frac{\pi}{2} & \Psi = 0 & \dfrac{\partial \Psi }{ \partial y} = 0 \\ & y = - \frac{\pi}{2} & \Psi = 0 & \dfrac{\partial \Psi }{ \partial y} = 0 \end{array}$$

where $$\Delta^2$$ is the biharmonic operator.

$$\Delta ^2 \Psi = \frac{\partial ^4 \Psi }{\partial x^4} + 2 \frac{\partial^4 \Psi }{\partial x^2 \partial y^2} + \frac{\partial ^4 \Psi }{\partial y^4}$$

I looked into DEigenvalues but I couldn't understand how the BCs are determined.

Can anyone help me to tell Mathematica to find the first five eigenvalues of this problem for me? I am interested to use finite elements.

• mathematica.stackexchange.com/questions/135893/… Jun 21 '20 at 14:55
• @user21: Can I ask that why you closed this question? This is an eigen-value problem while in the linked question you addressed to solve an inhomogeneous biharmonic BVP. Jun 22 '20 at 8:03
• My understanding is that your question is the same as the one referenced. Why don't you update your question with an explanation of why the linked answer is not applicable to your question? If that is the case then we can reopen your question, I'd say. Jun 22 '20 at 8:08
• I see there might be some issue in applying the result from the linked question to NDEigenstem but it would be good to point out that that issue actually is? Jun 22 '20 at 8:16
• @xzczd: If you have some innovative approach in mind, you are welcomed to present it here. :) Jun 22 '20 at 14:14

If you have some innovative approach in mind…

Hmm… I have a solution based on finite difference method (FDM), which is said to be the earliest numeric method for PDE solving. Anyway, it might be viewed as innovative in the sence that it shows FEM isn't the only possible numeric approach, so let me post it here.

## FDM Based Solution

I've used pdetoae for the generation of difference equations.

lap = Laplacian[#, {x, y}] &;
domain = {-Pi/2, Pi/2};
togeneralpattern = (a_ -> b_) :> Quiet@((a /. v : x | y -> v_) -> b);
With[{u = u[x, y]}, lhs = lap@lap@u;]

points = 400;
grid = Array[# &, points, domain];
difforder = 4;
ptoafunc = pdetoae[u[x, y], {grid, grid}, difforder];
discretelhs = ptoafunc@lhs;

ptoafuncmid = pdetoae[mid[x], grid, difforder];
eqmid = mid'[#] == 0 & /@ domain // ptoafuncmid
varmid = mid /@ grid[[{2, -2}]]
rulebc1 = Flatten@Table[rulemid, {mid, {u[x, #] &, u[#, y] &}}] /. togeneralpattern

del = #[[3 ;; -3, 3 ;; -3]] &;
innerlhs = Flatten@del@discretelhs;
innervar = Flatten@del@Table[u[x, y], {x, grid}, {y, grid}];
Block[{u}, Set @@@ rulebc1;
{barray, marray} = CoefficientArrays[innerlhs, innervar]]; // AbsoluteTiming
(* {38.2806, Null} *)

{val, vec} = Eigensystem[marray // N, -10]; // AbsoluteTiming
(* {21.8352, Null} *)

mat = Partition[#, points - 4] & /@ vec;
val
(* {454.983, 454.983, 279.493, 279.493, 179.43,
177.74,  120.223, 55.2993, 55.2993, 13.2938} *)

ArrayPlot[#, ColorFunction -> "AvocadoColors"] & /@ mat


You may further adjust points and difforder to check if the result is reliable.

Remark

The Dirichlet b.c.s are essentially set by del, to be precise, by removing outermost elements of Table[u[x, y], {x, grid}, {y, grid}].

## Update: Pseudospectral Method Based Solution

If you've read some of my answers using pdetoae, you may have noticed I've directly substitute the Neumann b.c.s into the system with the Set @@@ rulebc1; here, which is an approach I seldom took to impose b.c., because it makes the coding more involved. Nevertheless, if I replace part of discretized PDEs with b.c.s as I've usually done, we'll have to calculate generalized eigenvalues AFAIK, and the underlying algorithm of Eigensystem for generalized eigenvalues seems to be much less efficient.

It's better to illustrate with an example. In the following solution, pseudospectral method, for which direct substitution of b.c.s doesn't work well because the expressions of discretized system are so large in this case, is implemented.

lap = Laplacian[#, {x, y}] &;
domain = {-Pi/2, Pi/2};
With[{u = u[x, y]}, eq = {lap@lap@u, u};
{bcx, bcy} = Table[{{u, ϵ u}, {D[u, var], ϵ u}} /.
var -> boundary, {var, {x, y}}, {boundary, domain}]];

Off[NDSolveFiniteDifferenceDerivative::ordred]
points = 35;
CGLGrid[{xl_, xr_}, n_] := 1/2 (xl + xr + (xl - xr) Cos[(π Range[0, n - 1])/(n - 1)]);
grid = N@CGLGrid[domain, points];
difforder = points;
ptoafunc = pdetoae[u[x, y], {grid, grid}, difforder];
del = #[[3 ;; -3]] &;
discreteeq = del /@ del@# & /@ ptoafunc@eq;
discretebcx = Map[del, ptoafunc@bcx, {3}];
discretebcy = ptoafunc@bcy;
discretesys =
Flatten@{discreteeq[[#]], discretebcx[[All, All, #]], discretebcy[[All, All, #]]} &;
var = Table[u[x, y], {x, grid}, {y, grid}] // Flatten;
{barray, marray} = CoefficientArrays[discretesys@1, var]; // AbsoluteTiming
(* {2.33509, Null} *)
{barray2, aarray} = CoefficientArrays[discretesys@2 /. ϵ -> 0., var]; // AbsoluteTiming
(* {0.0019331, Null} *)
{val, vec} = Eigensystem[{marray, aarray} // N, -10]; // AbsoluteTiming
(* {11.0801, Null} *)
MaxMemoryUsed[]/1024^3. GB
(* 0.581552 GB *)
val
(*
{454.983 + 0. I, 454.983 + 0. I, 279.493 + 0. I, 279.493 + 0. I, 179.43 + 0. I,
177.74 + 0. I,  120.223 + 0. I, 55.2993 + 0. I, 55.2993 + 0. I, 13.2938 + 0. I}
*)


As we can see, with only 35 grid points, we obtain a solution that's almost the same as that obtained with 4th order difference scheme and 400 grid points, which shows the power of pseudospectral method. However, just increase points a little to, say, 50, you'll see the timing of Eigensystem increase to 114.809 and memory usage to 2.61183 GB, while with the method in first section, the timing is 0.13103 s and memory usage 0.153116 GB. Changing difforder to 4 doesn't help. The efficiency of generalized eigenvalue calculation is a bottleneck here.

BTW, to visualize the eigenfunction, we can't simply use ArrayPlot in this case, because the grid is no longer uniform. Interpolation is necessary:

mat = Partition[#, points] & /@ vec;
func = ListInterpolation[#, {grid, grid}] & /@ mat;

(f \[Function]
DensityPlot[f[x, y], {x, ##}, {y, ##}, PlotPoints -> points,
ColorFunction -> "AvocadoColors"] & @@ domain) /@ func

• (+1), Thanks for the nice and comprehensive answer. :) Jun 23 '20 at 5:03
• @HoseinRahnama There's a mistake in my original implementation and now it's corrected. It's somewhat surprising the mistake doesn't influence the result that much, which seems to be consistent with the observation here. Jun 23 '20 at 5:31
• @xzczd why not use the spectral method? Jun 23 '20 at 20:41
• @ABCDEMMM If you mean Pseudospectral, with the technique mentioned in the comments here it can be used, but I'm afraid it's not the best choice in this case. I've added an implementation of pseudospectral method to my answer , you may play with it. Jun 24 '20 at 9:00
• @xzczd 我想知道回答中的图形是哪个具体的调和方程的数值解，能否给出具体的方程和边界条件。 Jun 24 '20 at 9:06

eqn = {Laplacian[u[x, y], {x, y}] - v[x, y],
Laplacian[v[x, y], {x, y}]};
bcs = DirichletCondition[u[x, y] == 0, True];
ufun = NDEigenvalues[{eqn, bcs}, {u, v}, {x, -Pi/2, Pi/2}, {y, -Pi/2,
Pi/2}, 5]


This will give a warning message about a possible non positive definite matrix. Which I think is OK, but I have not thought about this deeply.

{-3.05309*10^-15, -0.999828, -0.999979, -1.99933, -2.00264}


One other thing to point put is that I used DirichletCondition in u and not the NeumannValue on u. In FEM you can have the one or the other but not both at the same position. Note, however, that you could apply a different bc in v in this case, which uses an implicit Neumann zero.

• I'm sorry, but the result is undesired. The Neumann b.c. chosen in this case is $\frac{\partial^3 \Psi }{\partial x^3} = 0$. And doesn't seem to be an easy way to implement OP's b.c. with FiniteElement of Mathmatica at the moment. Jun 22 '20 at 10:17
• @xzczd and this is what I said in my text - did I not? Jun 22 '20 at 10:36
• Oh… sorry, I should have read the commentary more carefully. Jun 22 '20 at 10:51
• @user21 is "pdetoae" build-in subroutine? Jun 23 '20 at 20:54
• @ABCDEMMM user21 does suggest me to move the function to wolfram resource function repository here, but as mentioned there, the documenting of pdetoode/pdetoae` is non-trivial. Jun 24 '20 at 1:49