# Solution of PDE system does no satisfy boundary conditions

I have a problem with solution of the 2D PDE system. It appears that solution does not match boundary condition at all.

I've tried using separately generated mesh, and it had not helped, just making error oscillation frequency larger, but not lowering the error amplitude.

Adding "Method -> {"ExplicitRungeKutta", "SpatialDiscretization" -> {"FiniteElement", "MeshOptions" -> {"MaxCellMeasure" -> 0.00001}}}" does not make any difference.

The region is a rectangle (I tried oval too) with two round holes. There is no such problem in the case of just one hole.

Do you have any ideas making it better?

    system = {nxU[x, y] == D[HU[x, y], x] + (D[nyU[x, y], x, y] + D[nxU[x, y], x, x]),
nyU[x, y] == D[HU[x, y], y] + (D[nyU[x, y], y, y] + D[nxU[x, y], x, y]),
D[nyU[x, y], y] + D[nxU[x, y], x] == D[HU[x, y], y, y] + D[HU[x, y],x,x]};

FuncLst = {nxU[x, y], nyU[x, y], HU[x, y]};

n0=0.3;hIn0=-0.5;rOut=30;r0=1.;Centre1={-10,0};Centre2={10,0};

RegionOut=Rectangle[{-rOut,-rOut/2},{rOut,rOut/2}];
BoundaryOut=RegionBoundary[RegionOut];

RegionIn1=Disk[Centre1,r0];
RegionIn2=Disk[Centre2,r0];
RegionIn=RegionUnion[RegionIn1,RegionIn2];
BoundaryIn1=RegionBoundary[RegionIn1];
BoundaryIn2=RegionBoundary[RegionIn2];
BoundaryIn=RegionBoundary[RegionIn];

RegionAll=RegionUnion[RegionDifference[RegionOut,RegionIn],BoundaryIn];

n0func1x[x_,y_]=n0 (x-Centre1[[1]])/r0;
n0func1y[x_,y_]=n0 (y-Centre1[[2]])/r0;
n0func2x[x_,y_]=n0 (x-Centre2[[1]])/r0;
n0func2y[x_,y_]=n0 (y-Centre2[[2]])/r0;

BndCond={DirichletCondition[nxU[x,y]==0,{x,y}\[Element]BoundaryOut],
DirichletCondition[nyU[x,y]==0,{x,y}\[Element]BoundaryOut],
DirichletCondition[HU[x,y]==0,{x,y}\[Element]BoundaryOut],
DirichletCondition[nxU[x,y]==n0func1x[x,y],{x,y}\[Element]BoundaryIn1],
DirichletCondition[nyU[x,y]==n0func1y[x,y],{x,y}\[Element]BoundaryIn1],
DirichletCondition[nxU[x,y]==n0func2x[x,y],{x,y}\[Element]BoundaryIn2],
DirichletCondition[nyU[x,y]==n0func2y[x,y],{x,y}\[Element]BoundaryIn2],
DirichletCondition[HU[x,y]==hIn,{x,y}\[Element]BoundaryIn1],
DirichletCondition[HU[x,y]==hIn,{x,y}\[Element]BoundaryIn2]};

solPDE = NDSolve[Join[system, (BndCond /. hIn -> hIn0)], Join[FuncLst],{x, y} \[Element] RegionAll][[1]];
nxUs[x_, y_] = Evaluate[nxU[x, y] /. solPDE];
xs[t_] = Centre1[[1]] + Cos[t];
ys[t_] = Centre1[[2]] + Sin[t];
Plot[{nxUs[xs[t], ys[t]], n0func1x[xs[t], ys[t]]}, {t, 0, 2 \[Pi]}]


• Are you sure that you cannot simplify your example a bit to showcase your problem? Commented May 24, 2018 at 18:33
• You have a minor error in your code: in the fifth line from the bottom, you have FunLst when I think you mean FuncLst. Commented May 24, 2018 at 20:46
• @MichaelSeifert , thank you! I've edited the code. Commented May 24, 2018 at 20:49
• @MarcoB I understand that it's too large. But unfortunately I can not see thу way simplify it. Thу main code is on the solution region definition and seting the soundary condition. Commented May 24, 2018 at 20:51
• Which version of Mathematica are you using? Have you seen this question & answer? Commented May 24, 2018 at 21:07

This is not an answer (yet), just a bunch of observations and suspicions that will hopefully lead a savvier person than myself to find a solution.

The problem, as far as I can tell, stems from the InterpolationOrder of the results returned by NDSolve. By default, the solutions returned by NDSolve seem to be second-order interpolating functions in both directions:

Needs["DifferentialEquationsInterpolatingFunctionAnatomy"];

(* {2, 2} *)


What's odd is that NDSolve is supposed to allow you to specify the InterpolationOrder of the function it returns; but for this problem (maybe for all finite-element methods?), NDSolve just seems to ignore this option:

solPDE2 = NDSolve[Join[system, (BndCond /. hIn -> hIn0)], Join[FuncLst], {x, y} \[Element] mesh, InterpolationOrder -> 1][[1]];

(* {2, 2} *)


OK... let's force a first-order mesh instead, which should yield a first-order interpolating function:

Needs["NDSolveFEM"]
mesh = ToElementMesh[RegionAll, "MeshOrder" -> 1];
solPDE3 = NDSolve[Join[system, (BndCond /. hIn -> hIn0)],
Join[FuncLst], {x, y} \[Element] mesh][[1]];

(* {1, 1} *)


Excellent. What do the boundaries look like?

nxUsnew[x_, y_] = Evaluate[nxU[x, y] /. solPDE3];
Plot[{nxUsnew[xs[t], ys[t]], n0func1x[xs[t], ys[t]]}, {t, 0, 2 \[Pi]}]


.... OK then.

We can see that at most of the actual points on the mesh, the boundary conditions are pretty well satisfied. This is probably true for the original function as well: the result of NDSolve is in contact with the true boundary function at several points, presumably the mesh points; and it's only the regions between the mesh points (where interpolation is necessary) that disagree. Also, running this version of the code with "MeshOrder" -> 2 yields something similar (though not identical) to your original result.

However, I have no idea why the points at $\theta = 0$ and $\theta = \pi$ (i.e., when $y = 0$) are off of the curve. On my system (Mac OS, MM 10.4) these points are aberrant both when using the ToElementMesh method and when using your original method to find the solution. However, your original graphs don't seem to have this problem. Perhaps you might have better luck using the ToElementMesh method on your system.

Still, it would be helpful to know the answers to two questions:

• Is there a way to get NDSolve to use a first-order mesh on a finite-element problem, without resorting to loading additional packages?
• What's up with the points at $y = 0$ being off the curve, regardless of the method used?
• Thank you! That's still mistery to me, but setting "MeshOrder" -> 1 worked for me (Win, Ver. 11.1), while "MeshOrder" -> 2 leads to similar, even worse errors on the boundary (though default value is 2 and ToElementMesh and NDSolve should produce the same result). Unfortunately "MeshOrder" -> 1 makes difficult using solution derivatives. Commented May 24, 2018 at 23:11