# 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:

$$\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 May 30 '17 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? – user21 May 30 '17 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 May 30 '17 at 13:52
• Do you have a model problem to which you happen to have the solution? – user21 May 30 '17 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}] – xzczd Jun 6 '17 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,
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][];
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.

• Why the difference in results between FEM and FDM? try this please Plot[sol[0, y], {y, 0, ymax}] for both. – zhk Jun 3 '17 at 8:53
• @zhk ……I think user21 has made a mistake when deducing $c$. I've commented below his answer. – xzczd Jun 3 '17 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 Jun 3 '17 at 9:25
• @zhk Did you notice the right hand side of the equation in the paper is -1`? – xzczd Jun 3 '17 at 9:31
• Oh yes but the velocity potential cannot be negative overhere. I think something wrong with my PDE. – zhk Jun 3 '17 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. – xzczd Jul 10 '20 at 3:53