# Solving a Modified Biharmonic Equation on a Square

I am seeking to solve the differential equation

$$$$\left[\partial_{\overline{x}}^{4}+2\left(1+\delta\right)\partial_{\overline{x}}^{2}\partial_{\overline{y}}^{2}+\partial_{\overline{y}}^{4}\right]\psi=w^{2}\gamma\left(\partial_{\overline{x}}^{2}+\partial_{\overline{y}}^{2}\right)\psi$$$$

on a square geometry $$(\overline{x},\overline{y})\in[-1/2,1/2]^2$$. The boundary conditions on the square edges are expressed in the minimal (non-)working example of Mathematica code below:

reg = Rectangle[{-.5, -.5}, {.5, .5}];

dd[x_, a_] := UnitBox[x/a];
f[x_, a_, d_] = dd[x - d, a] - dd[x + d, a];

bi2[w_, γ_, δ_] := {
Derivative[4, 0][ψ][x, y] +
2*(1 + δ)*Derivative[2, 2][ψ][x, y] +
Derivative[0, 4][ψ][x, y] - w^2 * γ * (Derivative[2, 0][ψ][x, y] + Derivative[0, 2][ψ][x, y])
};

BCs[a_, d_] := {
Derivative[1, 0][ψ][x, -.5] == f[x, a, d],
Derivative[1, 0][ψ][x, .5] == f[x, a, d],
Derivative[0, 1][ψ][-.5, y] == f[y, a, d],
Derivative[0, 1][ψ][.5, y] == f[y, a, d]
};

ψSol = ParametricNDSolveValue[
{
bi2[w2, γ2, δ2] == 0,
BCs[a2, d2]
},
ψ, {x, y} ∈ reg, {a2, d2, w2, γ2, δ2},
Method -> {"PDEDiscretization" -> {"MethodOfLines",
"SpatialDiscretization" -> "TensorProductGrid"}}
]

ψSol[.01, .1, 1., .4, .2][.1, .1]


If I try NDSolve with the method set to "FiniteElement", Mathematica complains about the degree of the PDE being greater than 2, but if I use another method (as in the above code), Mathematica complains about the boundary conditions: "Boundary values may only be specified for one independent variable."

I would greatly appreciate any advice about how to tackle this problem with Mathematica. If it helps, I don't need the full solution in the square, just values of second derivatives at the center of the square, e.g. $$\psi^{\left(2,0\right)}\left(0,0\right)$$, $$\psi^{\left(1,1\right)}\left(0,0\right)$$.

EDIT:

I've attempted a slightly different approach by rewriting the differential equation as $$\left(\nabla^{2}\right)^{2}\psi-w^{2}\gamma\nabla^{2}\psi+2\delta\partial_{\overline{x}}^{2}\partial_{\overline{y}}^{2}\psi=0$$ and following the advice given in the following comment: https://mathematica.stackexchange.com/a/185530/53559. Namely, by introducing auxiliary functions, we can make the PDE "look" second order. We can also use the same trick to make the boundary conditions no longer contain derivatives. However, despite these tricks my difficulties remain unresolved. Here again is a minimal (non-)working example:

reg = Rectangle[{-.5, -.5}, {.5, .5}];

dd[x_, a_] := (Pi * a)^(-1) * 1/(1 + (x/a)^2);
f[x_, a_, d_] = -dd[x + d, a] + dd[x - d, a];

bi2[w_, γ_, δ_] := {
Derivative[2, 0][ψ][x, y] + Derivative[0, 2][ψ][x, y] - ψ2[x, y],
Derivative[1, 1][ψ][x, y] - ψ3[x, y],
Derivative[2, 0][ψ2][x, y] + Derivative[0, 2][ψ2][x, y]
+ 2 * δ * Derivative[1, 1][ψ3][x, y] - w^2 * γ * ψ2[x, y],
D[ψ[x, y], x] - ψx[x, y],
D[ψ[x, y], y] - ψy[x, y]
};

BCs[a_, d_] := {
ψx[x, -.5] == f[x, a, d],
ψx[x, .5] == f[x, a, d],
ψy[-.5, y] == f[y, a, d],
ψy[.5, y] == f[y, a, d]
};

ψSol = ParametricNDSolveValue[
{
bi2[w2, γ2, δ2] == {0, 0, 0, 0, 0},
BCs[a2, d2]
},
{ψ, ψ2, ψ3, ψx, ψy}, {x, y} ∈ reg, {a2, d2, w2, γ2, δ2}
];

ψSol[.01, .1, 1., .4, .2][[1]][.2, .1]


Note that I've also changed the boundary condition functions from box functions to Lorentzians, in response to comments to my original question.

The above code runs, and I get no complaints about the order of the PDE or the boundary conditions, but now the returned function $$\psi$$ is identically zero over the square. I feel like this is progress toward a working solution, but I would greatly appreciate any insight on how this might be modified to work.

• Plot[f[x, -1/2, 1/2], {x, -1, 1}] clearly shows the discontinuities at $x= \frac 1 4,\,x=-\frac 1 4$. Does that problem make sense at all? – user64494 Oct 11 '20 at 5:59
• Yes, these are fairly standard boundary conditions as far as I'm aware. But in the limit $a\to 0$, all the really matters is that the function dd[x,a] is a nascent delta function, i.e. the smoother Lorentzian dd[x_, a_] := (Pi*a)^(-1)*1/(1 + (x/a)^2) would be equivalent for my purposes if that helps. – sferics Oct 11 '20 at 11:00
• @spherics: You wrote "Yes, these are fairly standard boundary conditions as far as I'm aware". Can you give a reference to ground your statement? TIA. – user64494 Oct 11 '20 at 11:22
• I find your phrasing both inaccurate and unneccesarily hostile. If you are unfamiliar with box functions and their ubiquity in physics as boundary conditions, it strikes me as unlikely that you're well-equipped to address the question posed, rendering my continued indulgence of your pettiness without purpose. I've explained to you what box functions are and why they're used in this context - now please either address the question posed or refrain from commenting further. – sferics Oct 11 '20 at 15:38
• The problem is underdetermined. As shown in the post you've found, 8 b.c.s are needed to determine a particular solution of biharmonic equation. Your second sample gives a solution because the default setting of FiniteElement method is zero Neumann value, to be more specific, NDSolve has silently used zero Neumann value for ψ at the boundary in this case. – xzczd Oct 12 '20 at 3:27