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?

Improve this question
$\endgroup$
1
  • $\begingroup$ Please provide definition of mesh too. $\endgroup$ Commented Feb 25, 2024 at 21:29

1 Answer 1

Reset to default
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.

Improve this answer
$\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$ Commented Feb 27, 2024 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$ Commented Feb 28, 2024 at 9:54

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.