For example, after solving the Poisson equation $-\nabla^2 T = 1$ on the unit circle using NDSolveValue and the following code

eqn = -Laplacian[u[x, y], {x, y}] == 1;
sol = NDSolveValue[{eqn, DirichletCondition[u[x, y] == 0, True]}, 
  u, {x, y} \[Element] Disk[]]

how to

  • sample the solution at an arbitrary point in the domain.
  • integrate the "numerically obtained" temperature value over the region? This is particular relevant when solving a transient problem and one wants to know the amount of heat in the system.
  • integrate the heat flux over the boundary? In other words, find the gradient normal to the local boundary and do a line integration (surface integration in 3D) over the boundary. Mathematically, the flux is $F= -\int_{\partial \Omega } \nabla T\cdot \boldsymbol {n} \,dl$, where $\boldsymbol {n}$ is the local normal direction pointing outward of the domain.

I searched for a long time without finding a solution. Thanks!

  • $\begingroup$ To evaluate the function give it values for x and y, e.g., sol[1/2, 3/4] or sol[x, y] /. {x -> 1/2, y -> 3/4} To integrate over the Disk use NIntegrate[sol[x, y], {x, y} \[Element] Disk[]] I'll leave the rest for others. $\endgroup$
    – Bob Hanlon
    Commented Feb 18, 2023 at 22:55
  • $\begingroup$ @BobHanlon When trying the integration, MMA complained "Numerical integration converging too slowly; ..." Is there a way to control the accuracy? To me, the third question is significantly more complicated because one has to derive the gradient and identify the location and normal of a region's boundary. If you know how to do that, please post the solution as an answer and I will be happy to accept it. Thanks! $\endgroup$
    – Taozi
    Commented Feb 19, 2023 at 3:08
  • $\begingroup$ With v13.2.1 for Mac OS X ARM (64-bit), NIntegrate[sol[x, y], {x, y} \[Element] Disk[]] evaluates without any warning or error messages. $\endgroup$
    – Bob Hanlon
    Commented Feb 19, 2023 at 3:31
  • $\begingroup$ Hi, I am the developer for this functionality. I am curious to learn how you searched for this and did not find what you are looking for. Depending on what you say, I'd like to improve the find-ability of the information you found. $\endgroup$
    – user21
    Commented Feb 20, 2023 at 7:18
  • $\begingroup$ For your last question, there is an example about Joule heating that shows how to get various fluxes. You might also find the HeatTransfer monograph interesting. $\endgroup$
    – user21
    Commented Feb 20, 2023 at 7:22

1 Answer 1


i) Evaluate an interpolating function:

sol[0, 0]
(* 0.25 *)

The first example of the NDSolve ref page has an example of a variant of this; many more examples are in the InterpolatingFunction ref page.

ii) Integration over the region:

mesh = sol["ElementMesh"];
NIntegrate[sol[x, y], {x, y} ∈ mesh]

(* 0.392702 *) 

iii) You could use this:

 Grad[sol[x, y], {x, y}] . 
  BoundaryUnitNormal[x, y], {x, y} ∈ ToBoundaryMesh[mesh]]
(* -3.13394 *)

If you want to integrate only over a part of the boundary, just create a boundary mesh for that part. See the ToBoundaryMesh ref page or the ElementMesh generation tutorial. One last note, you can use things like

 NeumannValue[1, x >= 1/2], {x, y} ∈ ToBoundaryMesh[mesh]]

Let me give an example. We use this from the documentation as a starting point. There is enough text that explains the model and I am just copying it here. The interesting part will be below.

vars = {T[x, y], {x, y}};
Ω = Rectangle[{0, 0}, {0.02, 0.01}];
pars = <|"ThermalConductivity" -> 3|>;
BCTemp = 
  HeatTemperatureCondition[x == 0 || x == 0.02, vars, 
   pars, <|"SurfaceTemperature" -> 1173|>];
BCconvective = 
  HeatTransferValue[y == 0.01, vars, 
   pars, <|"AmbientTemperature" -> 323, 
    "HeatTransferCoefficient" -> 50|>];
BCradiation = 
  HeatRadiationValue[y == 0.01, vars, 
   pars, <|"AmbientTemperature" -> 323.|>];
eqn = {HeatTransferPDEComponent[vars, pars] == 
    BCconvective + BCradiation, BCTemp};
Tfun = NDSolveValue[eqn, T, {x, y} ∈ Ω];

For a plot see the documentation.

Now, we can use the boundary condition specification to compute the integrals over those parts of the boundary:

bmesh = ToBoundaryMesh[Tfun["ElementMesh"]]

NIntegrate[BCconvective /. T -> Tfun, {x, y} ∈ bmesh]
(* -681.075 *)

NIntegrate[BCradiation /. T -> Tfun, {x, y} ∈ bmesh]
(*-1162.59 *)

I am not 100% sure, but I think this is what you are looking for. Your flux F is exactly the NeumannValue.

  • $\begingroup$ Wondferful! The function BoundaryUnitNormal perhaps deserves better exposure. One last question, could you please explain your last example? I only saw the use of NeumannValue as boundary condition specification. Thanks. $\endgroup$
    – Taozi
    Commented Feb 21, 2023 at 23:23
  • $\begingroup$ Yes, BoundaryUnitNormal is new and not yet documented. Currently some boundary conditions can generate it (See the heat transfer BCs), If you have and can share an example where the BoundaryUnitNormal would be useful to you I can make the documentation happen. Also see this. Let me know what you think. $\endgroup$
    – user21
    Commented Feb 22, 2023 at 10:34
  • $\begingroup$ Finding the mass/momentum/heat flux is the main use I can think of. That's often what the engineers care the most about because it gives the drag force, heat/mass transfer efficiency, etc. A more mathematical use is to demo the divergence theorem by finding a region's volume via surface integration (e.g., area = NIntegrate[{x, 0} . BoundaryUnitNormal[x, y], {x, y} \[Element] bmesh]) in 2D. $\endgroup$
    – Taozi
    Commented Feb 26, 2023 at 16:09

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.