2
$\begingroup$

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?

$\endgroup$
1
  • $\begingroup$ Please provide definition of mesh too. $\endgroup$ Feb 25 at 21:29

1 Answer 1

1
$\begingroup$

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} =NDSolve`ProcessEquations[{op == 0 , DirichletCondition[f[x, y] == 1, x == xMax] }, 
f, {x, y} \[Element] Rectangle[]];

NDSolve`FEM`GetInactivePDE[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.

$\endgroup$
2
  • $\begingroup$ 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. $\endgroup$
    – Andrew L
    Feb 27 at 17:15
  • $\begingroup$ @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! $\endgroup$ Feb 28 at 9:54

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.