# NDSolve FEM-Solver fails for simple PDE

Trying to solve this pde (inspired by this question NeumannValue when differential equation does not use Laplacian )

pde = Derivative[2, 0][f][x, y] ==y^2*f[x, y] + NeumannValue[1, x == 1 ]
bc = DirichletCondition[f[x, y] == 1, x == 1]  ;
F = NDSolveValue[{pde, bc}, f, {x, 0, 1}, {y, 0, 1}] ;
Plot[{Derivative[1, 0][F][1, y] }, {y, 0, 1},GridLines -> {None, {1}}, PlotRange -> {0, 1},PlotLabel ->"NeumannValue[1,x\[Equal]1] isn't fullfilled!" ]


In contrast to that solution my "handmade Galerkin method" gives a smooth result, which fullfills the Neumann condition quite well

(If of interest here my code handmade Galerkin method )

I'm aware that NeumanValue in the first example is ignored because of DirichletCondition at the same boundary.

My question :

Is it possible to make FEM solution work for the case that Neumann and Dirichlet bcs are defined at the same place?

Thanks!

• Neumann and Dirichlet bcs are defined at the same place? from Neumann condition and Dirichlet condition at the same point one of the answers says applying Dirichlet boundary conditions will override your Neumann boundary conditions in the case of the finite element method Commented Mar 1 at 12:28
• @Nasser Thanks for your reply. Galerkin method gets by without this restriction I think. First partial integration "processes" the Neumann condition, then the unknowns of Dirichlet condition are inserted. This leads to the final set of equations. Commented Mar 1 at 13:09
• I would like to see the code. Also, note that Dirichlet conditions apply at points while NeumannValues apply at edges (in 1D these are also points) Commented Mar 1 at 16:51
• @user21 Thanks, I try to share the file but am unsure how to do it... Commented Mar 1 at 16:54
• @AlexTrounev Thank you for your reply. I only introduced Gaussian quadratur to make code faster. But it seems to work. Meanwhile I try(actual without success ;-) ) to use Mathametica DiscretizePDE functionality to make code run faster... Commented Mar 1 at 20:05

To verifier solution computed with Galerkin method we can use the Euler wavelets collocation method from my answer here. We add boundary conditions at y==0 and y==1 in a form of homogenous Neumann condition or in a form of pde. First version

UE[m_, t_] := EulerE[m, t];
psi[k_, n_, m_, t_] :=
Piecewise[{{2^(k/2) UE[m, 2^k t - 2 n + 1], (n - 1)/2^(k - 1) <= t <
n/2^(k - 1)}, {0, True}}];
PsiE[k_, M_, t_] :=
Flatten[Table[psi[k, n, m, t], {n, 1, 2^(k - 1)}, {m, 0, M - 1}]]
k0 = 3; M0 = 4; With[{k = k0, M = M0},
nn = Length[Flatten[Table[1, {n, 1, 2^(k - 1)}, {m, 0, M - 1}]]]];
dx = 1/(nn); xl = Table[l*dx, {l, 0, nn}]; zcol =
xcol = Table[(xl[[l - 1]] + xl[[l]])/2, {l, 2, nn + 1}]; Psijk =
With[{k = k0, M = M0}, PsiE[k, M, t1]]; Int1 =
With[{k = k0, M = M0}, Integrate[PsiE[k, M, t1], t1]];
Int2 = Integrate[Int1, t1];
Psi[y_] := Psijk /. t1 -> y; int1[y_] := Int1 /. t1 -> y;
int2[y_] := Int2 /. t1 -> y;
M = nn;

U1 = Array[a, {M, M}]; U2 = Array[b, {M, M}]; G1 =
Array[g1, {M}]; G2 = Array[g2, {M}]; G3 = Array[g3, {M}]; G4 =
Array[g4, {M}];

u1[x_, z_] := int2[x] . U1 . Psi[z] + x G1 . Psi[z] + G4 . Psi[z];
u2[x_, z_] := Psi[x] . U2 . int2[z] + z G2 . Psi[x] + G3 . Psi[x];
uz[x_, z_] := Psi[x] . U2 . int1[z] + G2 . Psi[x];
ux[x_, z_] := int1[x] . U1 . Psi[z] + G1 . Psi[z];
uxx[x_, z_] := Psi[x] . U1 . Psi[z];
uzz[x_, z_] := Psi[x] . U2 . Psi[z];
(*Equations and solution*)

eqn = Join[
Flatten[Table[(uxx[xcol[[i]], zcol[[j]]] -
zcol[[j]]^2 u1[xcol[[i]], zcol[[j]]]), {i, M}, {j, M}]],
Flatten[Table[
u1[xcol[[i]], zcol[[j]]] - u2[xcol[[i]], zcol[[j]]], {i, M}, {j,
M}]]]; bc =
Join[Flatten[Table[ux[1, zcol[[j]]] - 1 == 0, {j, M}]],
Flatten[Table[u1[1, zcol[[j]]] - 1 == 0, {j, 1, M}]],
Flatten[Table[uz[xcol[[j]], 0] == 0, {j, M}]],
Flatten[Table[uz[xcol[[j]], 1] == 0, {j, M}]]]; var =
Join[Flatten[U1], Flatten[U2], G1, G2, G3, G4];

eq = Join[Table[eqn[[i]] == 0, {i, Length[eqn]}],
bc];

{v, m} = CoefficientArrays[eq, var];

sol1 = LinearSolve[m // N, -v];

(*Visualization*)
rul = Table[var[[i]] -> sol1[[i]], {i, Length[var]}];

Plot3D[Evaluate[u1[x, y] /. rul], {x, 0, 1}, {y, 0, 1},
PlotLegends -> Automatic, ColorFunction -> "SunsetColors",
PlotRange -> All, Exclusions -> None]

Plot3D[Evaluate[u2[x, y] /. rul], {x, 0, 1}, {y, 0, 1},
PlotLegends -> Automatic, ColorFunction -> "SunsetColors",
PlotRange -> All, Exclusions -> None, PlotTheme -> "Marketing",
MeshStyle -> White]


In second version we use pde as follows

eq0 = Join[Table[uxx[xcol[[i]], 0] == 0, {i, M}],
Table[uxx[xcol[[i]], 1] - u1[xcol[[i]], 1] == 0, {i, M}]];

eq1 = eq =
Join[Table[eqn[[i]] == 0, {i, Length[eqn]}], Take[bc, 2 M],
eq0]; sol2 = CoefficientArrays[eq1, var]

sol = LeastSquares[sol2[[2]] // N, -sol2[[1]]]


Solution sol is not differ from shown above sol1, which is already obvious.

Update 1. We can compare solution u1[x,y]/.rul and u1[x,y]/.rul1 with

sol3 = NDSolveValue[{Derivative[2, 0][f][x, y] ==
y^2*f[x, y], {f[1, y] == 1, Derivative[1, 0][f][1, y] == 1}},
f, {x, 0, 1}, {y, 0, 1}]


We have

{Plot[Table[Abs[sol3[x, y] - u1[x, y] /. rul1], {x, 0, 1, .1}] //
Evaluate, {y, 0, 1}, FrameLabel -> {"y", "\[Delta]f"},
PlotLegends -> Automatic, Frame -> True,
PlotLabel -> "Absolute error for pde at y=0,1"],
Plot[Table[Abs[sol3[x, y] - u1[x, y] /. rul], {x, 0, 1, .1}] //
Evaluate, {y, 0, 1}, FrameLabel -> {"y", "\[Delta]f"},
PlotLegends -> Automatic, Frame -> True,
PlotLabel -> "Absolute error for NeumannValue[0,y==0||y==1]"]}


• This solution mifght be obtained with FEM-solver eqnsFEM = {Derivative[2, 0][f][x, y] == y^2*f[x, y] + NeumannValue[1, x == xMax(*&&0\[LessEqual]y\[LessEqual]1*)] + (1/2 E^-y (1 - E^(2 y) (-1 + y) + y)) NeumannValue[1, x == 0 ] , DirichletCondition[f[x, y] == 1, x == xMax] }; solFEM = NDSolveValue[eqnsFEM, f, {x, xMin, xMax}, {y, yMin, yMax}] ; too. As you can see there is no Neumanncondition at y==0||y==1 included. Commented Mar 3 at 17:01
• Interestingly the solution exactly fullfills NeumannValue and DirichletCondition. If we try to decrease the last NeumannCondition  .1 (1/2 E^-y (1 - E^(2 y) (-1 + y) + y)) NeumannValue[1, x == 0 ] this property is lost. Don't know why Commented Mar 3 at 17:08
• @UlrichNeumann Yes, you are right, FEM solution also computable, but homogeneous NeumannValue[0, y==0] is Automatic option in this case, since this is Automatic for any PDE without given boundary conditions. Commented Mar 3 at 18:04
• I'm not so sure about that because for this special pde only the derivative Derivative[1,0][f][x,y] at boundary x=const can be specified as NeumannValue. Commented Mar 3 at 21:24
• The solution also agrees with the analytical one F = Values@DSolve[{Derivative[2, 0][f][x, y] == y^2*f[x, y], f[1, y] == 1,Derivative[1, 0][f][1, y] == 1}, f, {x, 0, 1 }, {y, 0, 1}] [[1,1]], which only prescribes two boundary conditions. Commented Mar 4 at 8:20

This is not an answer but an extended comment. I think it would be good to simplify and speed up your code and make sure the FEM computations are correct. Following the FEM Programming tutorial we get this:

Needs["NDSolveFEM"]

pde = Derivative[2, 0][f][x, y] ==
y^2*f[x, y] + NeumannValue[1, x == 1];
bc = DirichletCondition[f[x, y] == 1, x == 1];
F = data[{pde, bc}, f, {x, 0, 1}, {y, 0, 1}];


Using

NDSolveValue @@ F


gives what you had. This makes sure that we use the same input for NDSolveValus and for the system matrix extraction. Next we extract the system matrices:

{dPDE, dBC, vd, sd, md} = NDSolveFEMProcessPDEEquations @@ F;
sm = dPDE["StiffnessMatrix"];


Next, we deploy the boundary conditions:

DeployBoundaryConditions[{lv, sm}, dBC];


Solve the linear equations:

result = LinearSolve[sm, lv];


And post process the solution:

NDSolveSetSolutionDataComponent[sd, "DependentVariables",
Flatten[result]];
ProcessPDESolutions[md, sd]


Now, this solve the equation in exactly the same way like NDSolve would solve the PDE. What we can do next, is ignore the DirichletCondition. This will give a warning but we want to extract the load vector as if the DirichletCondition did not overwrite it:

F = data[{pde(*,bc*)}, f, {x, 0, 1}, {y, 0, 1}];
{dPDE, dBC, vd, sd, md} = NDSolveFEMProcessPDEEquations @@ F;
smDummy = dPDE["StiffnessMatrix"];
DeployBoundaryConditions[{lv, smDummy}, dBC];
lv


Now, you have the load vector with the NeumannValue deployed. In your notebook you transform the PDE into an optimization problem and I have some difficulty to understand what you are doing there. If you could show this minimization approach based on the stiffness matrix and the NeumannValue, that would be useful, for me to understand how you want to solve this.

Now, let's do another experiment. Now we ignore the DirichletBoundary condition altogether:

pde = Inactive[Div][{{-1, 0}, {0, 0}} .
Inactive[Grad][f[x, y], {x, y}], {x, y}] - y^2*f[x, y] ==
NeumannValue[1, x == 1];
(*bc=DirichletCondition[f[x,y]==1,x==1];*)
F = data[{pde(*,bc*)}, f, {x, y} \[Element] mesh];


Proceed as before, now we get an expected message:

{dPDE, dBC, vd, sd, md} = NDSolveFEMProcessPDEEquations @@ F;
sm = dPDE["StiffnessMatrix"];


DeployBoundaryConditions[{lv, sm}, dBC];
result = LinearSolve[sm, lv];
NDSolveSetSolutionDataComponent[sd, "DependentVariables",
Take[Flatten[result], md["DegreesOfFreedom"]]];
{ifun} = ProcessPDESolutions[md, sd];


Now, we look at the result in the same way the OP has done:

Plot[Evaluate[{Derivative[1, 0][ifun][1, y],
Derivative[1, 0][ifun][0, y]}], {y, 0, 1}]


• Thanks for your effort. See my answer(extended comment) which hopefully clarifies my approach. Commented Mar 4 at 8:34
• Is data in the definition of F==data[...] a protected Mathematica variable? Commented Mar 4 at 8:39
• No, why should it be? See update. Commented Mar 4 at 9:34
• Because data is undefined in the code , I thought. Commented Mar 4 at 9:37
• Command ProcessPDESolutions[md, sd] in your code doesn't return an interpolation object on Mathematica v12.2. Commented Mar 4 at 9:40

modified

Meanwhile I could DiscretizePDE too. Here my code

{xMin, xMax, yMin, yMax} = {0, 1, 0, 1};
Needs["NDSolveFEM"]


meshing

mesh = ToElementMesh[Rectangle[{xMin, yMin}, {xMax, yMax}],
"MaxCellMeasure" -> 0.0001 ,
"MeshElementType" -> "TriangleElement", "MeshOrder" -> 1] ;
index = MeshCells[MeshRegion[mesh], 0][[All, 1]];

rand = mesh["BoundaryElements"][[All, 1]] // Flatten[#, 1] &;
indexR = Flatten[rand ] // DeleteDuplicates;
pi = mesh["Coordinates"]; (* alle Punkte *)
fi = Table[f[i], {i, 1, Length[pi]}];  (* Unbekannte DOFs *)

reg1 = ImplicitRegion[0 <= y <= yMax && x == xMax, {x, y}]; (* Randx\[Equal]xMax*)
rand1 = Select[rand ,RegionMember[reg1, pi[[ #[[1]]]]] &&RegionMember[reg1, pi[[ #[[2]]]]] &];(* an diesem Rand gilt f[xMax,y]\[Equal]1*)
indexR1 = Flatten[rand1] // DeleteDuplicates;


discretization

vd = NDSolveVariableData[{"DependentVariables","Space"} -> {{f}, {x, y}}];
sd = NDSolveSolutionData[{"Space"} -> {mesh}];
cdata = InitializePDECoefficients[vd, sd,
"DiffusionCoefficients" -> {{  {{-1, 0}, {0, 0}}}},
"ReactionCoefficients" -> {{-y^2} }];
bcdata = InitializeBoundaryConditions[vd,sd, {{NeumannValue[1, x == 1]}}   ]
mdata = InitializePDEMethodData[vd, sd];

(*Discretization*)
dpde = DiscretizePDE[cdata, mdata, sd] ;
{load, stiffness, damping, mass} = dpde["All"];

dbc = DiscretizeBoundaryConditions[bcdata, mdata, sd];


solution

RB1 = Map[fi[[#]] == 1 &, indexR1];
solfun = ElementMeshInterpolation[{mesh}, fi /. mini[[2]]];

Plot3D[ solfun[x, y] , {x, y} \[Element] mesh, AxesLabel -> {x, y} ,
PlotRange -> {{0, 1}, All}[[-1]], MeshFunctions -> (#3 &)]


Plot [ {Derivative[1, 0][solfun][1, y],Derivative[1, 0][solfun][0, y]} , {y, 0, 1}]


The result shows , the two Neumann-conditions and the DirichletCondition are fullfilled quite well.

It seems to be possible to get plausible results with Galerkin method, even if we take NeumannValue and Dirichletcondition at the same place of the boundary!

With the last example (@user21 's extended comment) the minimization problem with matrices smdummy and lv follows after introducing the dirichletcondition

RB1 = Map[fi[[#]] == 1 &, indexR1];(* all meshpoints x==1*)
mini = NMinimize[{(smdummy.fi + Flatten[lv]) . (smdummy.fi+    Flatten[lv]), RB1}, fi];

• @user21 But it gives the same result found in galerkinUN.nb. Commented Mar 4 at 9:46
• I'd say, because what you ask for can not be fulfilled. Commented Mar 4 at 10:13
• @user21 But the last plot shows the Neumannconditions are "roughly" fullfiled!? Commented Mar 4 at 10:14
• @user21 I am not yet convinced from your last comment, see my modified answer(extended comment) with refined mesh Commented Mar 4 at 12:15
• You can use this: bcdata2 = InitializeBoundaryConditions[vd, sd, {{DirichletCondition[f[x, y] == 1, x == 1]}}]; indexR1 = bcdata2["Stationary"][[1, 1, 3, 2, 1]] to compute the indexR1 Commented Mar 5 at 9:22