# How to input Robin boundary conditions for nonstandard Laplace equation?

What is the correct way to input these boundary conditions for the following nonstandard Laplace equation, whose coefficients of $$\frac{\partial^2 u}{\partial x^2}$$ and $$\frac{\partial^2 u}{\partial y^2}$$ term are not the same?

$$\epsilon^2\frac{\partial ^2u}{\partial x^2}+\frac{\partial ^2u}{\partial y^2}=1$$

$$\left.\frac{\partial u}{\partial x}\right|_{x=0}=0,\,\,\left.\frac{\partial u}{\partial y}\right|_{y=0}=0,\,\left.u=-2n\frac{\partial u}{\partial y}\right|_{y=1},\,\left.u=-2\epsilon n\frac{\partial u}{\partial x}\right|_{x=1}$$

Here is my try

ClearAll["Global*"]

Needs["NDSolveFEM"];
xmax = 1; ymax = 1;

epsilon = 0.5; n=0.1;

Ω := Rectangle[{0, 0}, {xmax, ymax}];

nv1 = NeumannValue[0, x == 0];
nv2 = NeumannValue[0, y == 0];
nv3 = NeumannValue[-n*epsilon*u[x, y], x == xmax];
nv4 = NeumannValue[-n*u[x, y], y == ymax];

sol = NDSolve[{epsilon^2*D[u[x, y], x, x] + D[u[x, y], y, y] == 1 + nv1 + nv2 + nv3 + nv4},
u, {x, y} ∈ Ω];


Here is something to get you started:

c = -{{eps^2, 0}, {0, 1}} /. eps -> epsilon;
op = Div[c.Grad[u[x, y], {x, y}], {x, y}];
sol = NDSolveValue[{op ==
1 - NeumannValue[u[x, y], x == xmax || y == ymax]},
u, {x, y} ∈ Ω]


You'd need to think about the signs. (Do you really want c to be positive). The point is you need to fit your equation into the coefficient form for PDEs for the finite element solver to work. Have a look at how the coefficients of the PDE are related to the coefficients in the NeumannValue. It's important to understand that the coefficients of the PDE are not independent of the coefficients of the NeumannValue. You can find more information in the section Partial Differential Equations and Boundary Conditions of the documentation. As an alternative the details section of the reference page for NeumannValue has additional information.

Update:

Let's assume your equation is:

$$\nabla\cdot (-c \nabla u) - 1=0$$

this implies that the Neumann/Robin operator is:

$$n \cdot (c \nabla u)=g + q u$$

In an initial update I miss read the boundary condition. Because there is an $\epsilon$ and not an $\epsilon^2$ in one of the boundary conditions we model the Robin condition by using NeumannValue one on each side.

a = 1/2; b = 1;
epsilon = a/b; n = 1/10;
Ω = Rectangle[{0, 0}, {xmax = 2 b, ymax = 2 a}];
c = -{{epsilon^2, 0}, {0, 1}};
op = Div[c.Grad[u[x, y], {x, y}], {x, y}] - 1;
g = 0; q = 1/(2 n);
solFEM = NDSolveValue[{op == -NeumannValue[g + epsilon*q*u[x, y],
x == xmax] - NeumannValue[g + q*u[x, y], y == ymax]},
u, {x, y} ∈ Ω];
Plot3D[solFEM[x, y], {x, 0, xmax}, {y, 0, ymax}]


This result agrees with the FDM version by @xzczd once PDE coefficients are made to match. Note that this approach, although it requires a bit of thought, needs much less code than the FDM version.

Comparison with FDM version:

With[{u = u[x, y]}, eq = epsilon^2 D[u, x, x] + D[u, y, y] == -1;
{bc@x, bc@y} = {{D[u, x] == 0 /. x -> 0,
u == -2 epsilon n D[u, x] /. x -> xmax}, {D[u, y] == 0 /. y -> 0,
u == -2 n D[u, y] /. y -> ymax}};]


Gives:

Plot3D[solFEM[x, y] - sol[x, y], {x, 0, xmax}, {y, 0, ymax}]


• First of all thank you. Is there otherway (NDSolve) to solve this problem without using FEM?
– zhk
Commented May 30, 2017 at 10:47
• Not directly, NDSolve will use the FEM for elliptic PDEs. You could, however, code your own finite difference with FiniteDifferenceDerivative; but I think it should be possible to solve this with FEM. What are your concerns? Commented May 30, 2017 at 13:16
• My main concern is the sign with c. I read the documentation and was unable to grasp it. Is it the same as multiplying the left hand of the PDE by a negative sign? or you are having c with a negative sign because of the robin boundary conditions?
– zhk
Commented May 30, 2017 at 13:52
• Do you have a model problem to which you happen to have the solution? Commented May 30, 2017 at 14:51
• Er…for example epsilon = 0.5; n = 1/10;c = {{coe@1, 0}, {0, coe@2}}; With[{u = u[x, y]}, eq@1 = epsilon^2 D[u, x, x] + D[u, y, y] == -1; eq@2 = Div[-c.Grad[u, {x, y}], {x, y}] == f; {bcx@1, bcx@2} = {u == -2 epsilon n D[u, x], {1, 0}.(c.Grad[u, {x, y}]) == q@1 u} /. x -> xmax; {bcy@1, bcy@2} = {u == -2 n D[u, y], {0, 1}.(c.Grad[u, {x, y}]) == q@2 u} /. y -> ymax;];FindInstance[Flatten[(Reap@Collect[Subtract @@ #@1 - k@# Subtract @@ #@2, _[x | xmax, y | ymax], Sow])[[-1,1]] & /@ {eq, bcx, bcy}] == 0 // Thread, Flatten@{{#@1, #@2} & /@ {coe, q}, k /@ {eq, bcx, bcy}, f}] Commented Jun 6, 2017 at 3:05

I'd like to add a solution based on finite difference method (FDM). The advantage of this solution is we don't need to manually transform the equation to any standard form. I'll use pdetoae for the generation of difference equation.

xmax = 1; ymax = 1;
epsilon = 0.5; n = 0.1;
With[{u = u[x, y]},
eq = epsilon^2 D[u, x, x] + D[u, y, y] == 1;
{bc@x, bc@y} = {{D[u, x] == 0 /. x -> 0, u == -2 epsilon n D[u, x] /. x -> xmax},
{D[u, y] == 0 /. y -> 0, u == -2 n D[u, y] /. y -> ymax}};]
domain@x = {0, xmax};
domain@y = {0, ymax};
points@x = 25;
points@y = 25;
difforder = 4;
(grid@# = Array[# &, points@#, domain@#]) & /@ {x, y};
var = Outer[u, grid@x, grid@y] // Flatten;

(* Definition of pdetoae isn't included in this post,
please find it in the link above. *)
ptoafunc = pdetoae[u[x, y], grid /@ {x, y}, difforder];
removeredundance = #[[2 ;; -2]] &;
ae = removeredundance /@ removeredundance@ptoafunc@eq;
aebc@x = removeredundance /@ ptoafunc@bc@x;
aebc@y = ptoafunc@bc@y;
solrule = Solve[{ae, aebc@x, aebc@y} // Flatten, var][[1]];
solpoints = N@solrule /. (u[x_, y_] -> value_) :> {x, y, value};
sol = Interpolation[solpoints]

(* The following is an alternative method for obtaining sol,
it's more challenging to understand, but more efficient. *)
(*
{b, m} = CoefficientArrays[{ae, aebc@x, aebc@y} // Flatten, var];
sollst = LinearSolve[m, -b];
sol = ListInterpolation[ArrayReshape[sollst, {points@x, points@y}], domain /@ {x, y}]
*)
Plot3D[sol[x, y], {x, 0, xmax}, {y, 0, ymax}]


# Update

If you still feel confused about removeredundance, the followings are 2 alternatives that don't require you to remove equations from the system:

fullsys = Flatten@ptoafunc@{eq, bc@x, bc@y};

(* Alternative 1: *)

lSSolve[obj_List, constr___, x_, opt : OptionsPattern[FindMinimum]] :=
FindMinimum[{1/2 obj^2 // Total, constr}, x, opt]
lSSolve[obj_, rest__] := lSSolve[{obj}, rest]

solrule = Last@
lSSolve[Subtract @@@ fullsys, var]; // AbsoluteTiming

solpoints = N@solrule /. (u[x_, y_] -> value_) :> {x, y, value};
sol = Interpolation[solpoints]

(* Alternative 2: *)

{b, m} = CoefficientArrays[fullsys, var];

sollst = LeastSquares[m, -b]; // AbsoluteTiming
sol = ListInterpolation[ArrayReshape[sollst, {points@x, points@y}], domain /@ {x, y}]


You can check this post to learn more about lSSolve.

# Update 2

With my allowfemdbc we can use FiniteElement method without manually setting up NeumannValue:

solfem = allowfemdbc@NDSolveValue[{eq, bc /@ {x, y}},
u, {x, 0, xmax}, {y, 0, ymax}];

Plot3D[solfem[x, y], {x, 0, xmax}, {y, 0, ymax}]


The solution looks the same as above, so I'd like to omit it here.

• Why the difference in results between FEM and FDM? try this please Plot[sol[0, y], {y, 0, ymax}] for both.
– zhk
Commented Jun 3, 2017 at 8:53
• @zhk ……I think user21 has made a mistake when deducing $c$. I've commented below his answer. Commented Jun 3, 2017 at 9:11
• But if we compare the FDM results with this paper dropbox.com/s/amvhn5xdq49lmir/… Both looks different. To get the same type of behaviour we need to have -sol[0,y].
– zhk
Commented Jun 3, 2017 at 9:25
• @zhk Did you notice the right hand side of the equation in the paper is -1`? Commented Jun 3, 2017 at 9:31
• Oh yes but the velocity potential cannot be negative overhere. I think something wrong with my PDE.
– zhk
Commented Jun 3, 2017 at 9:32

Version 12.1 can solve this exactly: Laplace equation with robin boundary conditions

• The method therein won't help. OP's equation is a non-standard Laplace equation whose coefficient is special. Please read user21's answer above carefully for more detail. Commented Jul 10, 2020 at 3:53