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?

  • $\begingroup$ 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? $\endgroup$
    – xzczd
    Commented Jun 12, 2019 at 5:50
  • $\begingroup$ 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. $\endgroup$
    – Alex T.
    Commented Jun 12, 2019 at 20:06
  • $\begingroup$ 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. ) $\endgroup$
    – xzczd
    Commented Jun 13, 2019 at 5:25

1 Answer 1


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

<< NDSolve`FEM`
\[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]\
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]



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.