# The context:

I am attempting to solve the differential equation d^2f(x,y)/dx^2=y^2*f(x,y) on the rectangle given by (x_min,x_max)=(0,1), (y_min,y_max)=(0,1) with the boundary conditions f(x_max,y)=1 and df/dx=-1 at (x_max,y). This can be done easily with a simple NDSolve call, e.g.

eqns = {D[f[x, y], x, x] == y^2*f[x, y],
f[xMax, y] == 1, (D[f[x, y], x] /. {x -> xMax}) == 1};
sol = NDSolve[eqns, f, {x, xMin, xMax}, {y, yMin, yMax}][[1]];


However, this is a toy model for a more complicated problem, which I need to solve on a predetermined coordinate grid mesh, e.g.,

xList = Range[xMin, xMax, 0.02];
yList = Range[yMin, yMax, 0.02];
coords = Flatten[Table[{x, y}, {x, xList}, {y, yList}], 1];
nPointsX = Length[xList];
nPointsY = Length[yList];
meshElements =
Flatten[Table[
QuadElement[{{(i - 1)*nPointsY + j, i*nPointsY + j,
i*nPointsY + j + 1, (i - 1)*nPointsY + j + 1}}], {i,
1, (nPointsX - 1)}, {j, 1, (nPointsY - 1)}], 1];
mesh = ToElementMesh["Coordinates" -> coords,
"MeshElements" -> meshElements];


(which creates a rectangular grid on the rectangular domain) as opposed to a grid that Mathematica generates internally. I attempt to do this by the following:

eqnsmesh = {D[f[x, y], x, x] ==
y^2*f[x, y] + NeumannValue[1, x == xMax],
DirichletCondition[f[x, y] == 1, x == xMax]};
solmesh = NDSolve[eqnsmesh, f, Element[{x, y}, mesh]][[1]];


but this produces a different solution for f[x,y]. Indeed, I found that attempting to adapt my original PDE to use the NeumannValue feature:

eqns2 = {D[f[x, y], x, x] ==
y^2*f[x, y] + NeumannValue[1, x == xMax],
DirichletCondition[f[x, y] == 1, x == xMax]};
sol2 = NDSolve[eqns2, f, {x, xMin, xMax}, {y, yMin, yMax}][[1]];


produces the same solution for f[x,y] as found with eqnsmesh and solmesh above.

# The question:

The Mathematica documentation for NeumannValue, while sparse, does indicate that this feature works on an understanding that the PDE uses the Laplacian, https://reference.wolfram.com/language/FEMDocumentation/tutorial/FiniteElementBestPractice.html#1529668360. My PDE, however, only uses two x derivatives, and therefore is not expressed in terms of the Laplacian. How can I format the NeumannValue arguments so that the solution produces the same function as eqns and sol, which don't require the NeumannValue feature to solve?

• Please provide definition of mesh too. Feb 25 at 21:29

modified

The NeumannValue for your problem is well defined.

{xMin, xMax , yMin, yMax} = {0, 1, 0, 1}; (* you didn't provide these values*)
op = Derivative[2, 0][f][x, y] == y^2*f[x, y];
{state} =NDSolveProcessEquations[{op == 0 , DirichletCondition[f[x, y] == 1, x == xMax] },
f, {x, y} \[Element] Rectangle[]];

NDSolveFEMGetInactivePDE[state]


$$\left\{\nabla _{\{x,y\}}\cdot \left(-\left( \begin{array}{cc} -1 & 0 \\ 0 & 0 \\ \end{array} \right).\nabla _{\{x,y\}}f(x,y)\right)-y^2 f(x,y)\right\}$$

Coefficient of Div .  (see FEM documentation), here

-{{-1,0},{0,0}}.Grad[f]=={Derivative[1,0]f][x,y],0}


, multiplied by the unitnormal {nx,ny} of the boundary

{Derivative[1,0]f][x,y],0}.{nx,ny}==Derivative[1,0]f][x,y]nx


gives the NeumannValue which Mathematica FEM-solver allows.

As you might see for your problem NeumannValue only at boundaries x==xMin and x==xMax might be prescribed.

• Thank you for the response! I'm still a bit confused on how NeumannValue would be implemented here. If I were to solve the PDE in the form you've given, I would use something like op==0+NeumannValue[ ? ]. What goes inside the NeumannValue function to specify the actual boundary conditions? From your answer, I see that in terms of the variables given in the documentation, c = {{-1,0},{0,0}}, a=-y^2, and alpha=beta=gamma=f=0. But the connection between c and the argument of the NeumannValue call that I need to implement is unclear, if such a connection exists. Feb 27 at 17:15
• @AndrewL From yanswer you get -{{-1,0},{0,0}}.Grad[f]=={Derivative[1,0]f][x,y],0}. This resulty has to be multiplied by the unit normal of the boundary to get the NeumanValue`. See my modified answer! Feb 28 at 9:54