# Variation of the biharmonic equation with Neumann conditions

I am currently writing a script to plot the solution of a variant of the biharmonic equation. In this case the equation I want to solve is

Laplacian[\[Alpha] Laplacian[u[x,y],{x,y}], {x,y}] + \[Beta]Laplacian[u[x,y], {x,y}] + \[Gamma] u[x,y] == 0


where \[Alpha], \[Beta], and \[Gamma] are prescribed constants. I have successfully been able to plot a solution using the following:

biharm = {Laplacian[u[x, y], {x, y}] == v[x, y],
Laplacian[\[Alpha] v[x, y], {x,y}] == -\[Beta] v[x, y] - \[Gamma] u[x, y]};
bound = {u[0, y] == 0, u[1, y] == 0, u[x, 0] == 0, u[x, 1] == 0,
v[0, y] == 0, v[1, y] == 0.5, v[x, 0] == 0, v[x, 1] == 0};
sol = NDSolveValue[{biharm, bound}, u, {x, 0, 1}, {y, 0, 1}];
Plot3D[sol[x, y], {x, 0, 1}, {y, 0, 1}, PlotRange -> All]


However, the boundary conditions with v[x,y] are not actually want I want to use; instead, I would like to prescribe conditions for D[u[x,y],x] and D[u[x,y],y] on the boundaries of my region as I have done in the first line of the boundary conditions.

When I try to implement this, I receive the following error:

The dependent variable in ... in the boundary condition DirichletCondition[...] needs to be linear.


From there, the script outputs a blank 3DPlot. The error seems to indicate that I would be able to specify something like NeumannValue instead of DirichletCondition, but none of the similar examples I have looked at seem to indicate how I can do that in my case. (For instance, I have seen NeumannValue[...] done with Laplace's equation, but I am unsure how to implement it here.)

Is there a way to do this here, or have I missed something in the other examples?

• 1. What boundary do you want to impose? 2. “When I try to implement this, I receive the following error…” Where's the corresponding code? – xzczd Jun 12 at 5:50
• Sorry, I want to impose this boundary instead: rather than v[0,y]==0, for instance, I would like (D[u[x, y], x] /. x -> 0)==0. (The remaining 3 are all changed in a similar fashion.) The corresponding code is exactly that as above, but with (e.g.) (D[u[x, y], x] /. x -> 0)==0 in place of those conditions on v[x,y]. I am not sure whether part of the issue is that I am providing 8 boundary conditions when I am providing two 2nd order PDEs, and maybe Mathematica expects four conditions each for u[x,y] and v[x,y]. I'm really not sure. Maybe there is a better way to input this 4th order PDE. – Alex T. Jun 12 at 20:06
• Then as you've noticed, you need to use NeumannValue, which should not be difficult for your case. (BTW here is a troublesome case. ) Just check the Details section of NeumannValue. (Another thing that's worth to mention is, zero NeumannValue can be omitted. ) – xzczd Jun 13 at 5:25

If you are only interested in the principle and not in solving a specific problem, then you can use such code.

<< NDSolveFEM
\[Alpha] = 1; \[Beta] = 1; \[Gamma] = 1;
reg = ImplicitRegion[0 <= x <= 1 && 0 <= y <= 1, {x, y}];
mesh = ToElementMesh[reg, MaxCellMeasure -> .001];
eq1 = -Laplacian[u[x, y], {x, y}] +
v[x, y]; eq2 = -Laplacian[\[Alpha] v[x, y], {x, y}] - \[Beta] v[x,
y] - \[Gamma] u[x, y];
bound = {DirichletCondition[v[x, y] == 0, y == 0],
DirichletCondition[u[x, y] == 0,
x == 0]};(*u[x,0]\[Equal]0,u[x,1]\[Equal]0,v[0,y]\[Equal]0,v[1,y]\
\[Equal]0.5,v[x,0]\[Equal]0,v[x,1]\[Equal]0};*)
sol = NDSolveValue[{eq1 ==
NeumannValue[1, x == 0] + NeumannValue[-1, y == 0],
eq2 == NeumannValue[0, True], bound}, u, {x, y} \[Element] mesh];

Plot3D[sol[x, y], {x, 0, 1}, {y, 0, 1}, PlotRange -> All,
AxesLabel -> Automatic, Mesh -> None, ColorFunction -> Hue]