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 Jun 12 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. Jun 12 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 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.

<< 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, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.