I am having issues in integrating the numerical solution obtained with FEM method and NDSolve, over a boundary. Specifically, i am integrating a specific combination of the gradient over a boundary obtaining errors which are much larger than the specified tolerances. I would like to have a way to avoid the errors and to have a faster integration.

Let me give you the full story:

I am solving the equation $\boldsymbol{\nabla}^2 f(x,y) -(0.001)^2 f(x,y)=0$ in two dimensions and in cartesian coordinates.

My domain consists in a very big outer circle of radius $L=50000$ and two small circles $C_1$ and $C_2$ of radius $R=1$ whose $x$ and $y$ center's coordinates are $\left(- \frac{dist}{2},0 \right)$ and $\left(\frac{dist}{2},0 \right)$, respectively. $dist$ is a parameter that I vary to fix the circles center to center distance.

On the external circle i have no flux boundary condition (Neumann boundary condition). Whereas on the two circles i have Dirchelet boundary condition:

$f(x,y)=1 +(x-\frac{dist}{2})^2-y^2$ on $C_2$

$f(x,y)=1 +(x+\frac{dist}{2})^2-y^2$ on $C_1$

This is the code that i use to generate the mesh and solve the problem (I specify the mesh elements on the boundaries):

L = 50000;
dist = 5.0;
nptspart = 2000;
nptsout = 60;
ptspart = 
 Table[{Cos[2.0*Pi*i/nptspart] - dist/2, Sin[2.0 Pi*i/nptspart]}, {i, 0, nptspart - 1}];
elempart = 
 Table[If[i < Length[ptspart], {i, i + 1}, {Length[ptspart], 1}], {i, 1, Length[ptspart]}];
ptspart1 = 
 Table[{Cos[2.0*Pi*i/nptspart] + dist/2, Sin[2.0 Pi*i/nptspart]}, {i, 0, nptspart - 1}];
elempart1 = 
 Table[If[i < Length[ptspart1], {nptspart + i, 
    nptspart + i + 1}, {nptspart + Length[ptspart1], 
    nptspart + 1}], {i, 1, Length[ptspart1]}];
ptsout = Table[{L Cos[2.0*Pi*i/nptsout], L Sin[2.0 Pi*i/nptsout]}, {i,
    0, nptsout - 1}];
elemout = 
 Table[If[i < Length[ptsout], {2 nptspart + i, 
    2 nptspart + i + 1}, {2 nptspart + Length[ptsout], 
    2 nptspart + 1}], {i, 1, Length[ptsout]}];
bmesh = ToBoundaryMesh[
  "Coordinates" -> Join[ptspart, ptspart1, ptsout], 
  "BoundaryElements" -> {LineElement[
     Join[elempart, elempart1, elemout]]}, "MeshOrder" -> 2];
mesh1 = ToElementMesh[bmesh, "MeshOrder" -> 2, 
  "RegionHoles" -> {{-dist/2, 0}, {dist/2, 0}}, 
  MeshQualityGoal -> "Maximal"];
ufun1 = NDSolveValue[{\!\(
\*SubsuperscriptBox[\(∇\), \({x, y}\), \(2\)]\(f[x, 
       y]\)\) - (0.001)^2 f[x, y] == 
    NeumannValue[0.0, x^2 + y^2 == L^2], 
    f[x, y] == 
     1 + ((x - dist/2)^2 - y^2), (x - dist/2)^2 + y^2 - 1^2 == 0], 
    f[x, y] == 
     1 + ((x + dist/2)^2 - y^2), (x + dist/2)^2 + y^2 - 1^2 == 0]}, 
  f, {x, y} ∈ mesh1, 
  Method -> {"PDEDiscretization" -> "FiniteElement"}];

The parameter nptspart defines how many points i have on the circles boundaries $C_1$ and $C_2$, whereas nptsout defines the number of points on the bigger circle of radius L.

I have now to perform an integral on the boundary of each circle:

$int_1=\int_{C_1} \left( \boldsymbol{\nabla} f \boldsymbol{\nabla} f -\frac{1}{2} \boldsymbol{\nabla} f \cdot \boldsymbol{\nabla} f \boldsymbol{I} \right) \cdot \boldsymbol{n} \; dl$

$int_2=\int_{C_2} \left( \boldsymbol{\nabla} f \boldsymbol{\nabla} f -\frac{1}{2} \boldsymbol{\nabla} f \cdot \boldsymbol{\nabla} f \boldsymbol{I} \right) \cdot \boldsymbol{n} \; dl$

Where $\boldsymbol{I}$ is the identity matrix in two dimensions, $\boldsymbol{n}$ the outwardly directed normal to the circles boundary, $\nabla f \nabla f$ is the dyadic product (also called Kronecker product), and $dl$ is the differential element on the circle line.

In orer to perform the integration, on the circle $C_2$ for example, i use the following code:

NIntegrate[(KroneckerProduct[Grad[ufun1[x, y], {x, y}], 
     Grad[ufun1[x, y], {x, y}]] - (1/
       2) (Grad[ufun1[x, y], {x, y}].Grad[ufun1[x, y], {x, y}])*
     IdentityMatrix[2]).{x - dist/2, y}, {x, y} ∈ 
  Circle[{dist/2, 0}, 1], AccuracyGoal -> 6, PrecisionGoal -> 6, 
 MaxRecursion -> 100, Method -> "InterpolationPointsSubdivision"]

However not only i get a lot of warning messages like:

NIntegrate::slwcon: Numerical integration converging too slowly; suspect one of the following: singularity, value of the integration is 0, highly oscillatory integrand, or WorkingPrecision too small.

NIntegrate::eincr: The global error of the strategy GlobalAdaptive has increased more than 400 times. The global error is expected to decrease monotonically after a number of integrand evaluations. Suspect one of the following: the working precision is insufficient for the specified precision goal; the integrand is highly oscillatory or it is not a (piecewise) smooth function; or the true value of the integral is 0. Increasing the value of the GlobalAdaptive option MaxErrorIncreases might lead to a convergent numerical integration. NIntegrate obtained 0.06632359368254395and 0.0005234490597411253 for the integral and error estimates.

but also the Integration is very slow, which is something I don't want as i will have to perform the integration, in my future research, for, say, 200 circles. therefore i need the integration to be fast on each boundary.

If i plot the results of the two integrals, for $dist=5$, as a function of the number of points for the circle boundaries (nptspart), I get something like this.

enter image description here

Where you can clearly see an error of the order of $\pm 0.001$ which is way bigger than the specified precision (6 digits) i specified in the integration.

Furthermore the error doesn't seem to reduce appreciably as the mesh on the circles gets refined.

POSSIBLE REASON: I think this has something to do with the interpolation function being evaluated in points which are not mesh points and are specified in the integration as:

  {x, y} ∈ Circle[{dist/2, 0}, 1]

Indeed if i let mathematica do a plot over the boundary $C_1$ (where the value of the function is specified by the Dirchelet boundary condition) i get weird results when i reduce the number of points on the boundary.

If i perform a parametric plot of the function $f(x(\theta),y(\theta))$ on $C_1$ (parametrised with the angle $\theta$ as shown in figure) for two different meshes

enter image description here

i get the following weird result:

enter image description here

In the above example i used directly the interpolation function obtained as result from the solution of NDSolve. An it is clear that is doing weird things on the boundary. So i guess that the gradient would be even worse.

So The question is. Is there a way to obtain the values of the function, and the gradient of the function, in the mesh nodes? Or if not, is it possible to perform an integration with the interpolation function that gives accurate results and is performed faster than it is at the moment? Do you have any tips?

  • $\begingroup$ your problem is not that the integral is done in a circle but the function evaluated on the elements. The problem is you are using fem for something is not made for. I guess the "real" problem you want to solve is two disks with a vanishing function in infinity. the result of the integral will not only depend on the nptspart but also the size of the system and your grid, making any parametric study impossible. The method to solve this is boundary element method, google bemlib for a good explanation. $\endgroup$ Commented Oct 24, 2016 at 17:43
  • 3
    $\begingroup$ I'm literally working on the same problem at this moment and have implemented BEM in mathematica for 2D in case you are interested on the code, I was going to post a question and answer it but still need to check the results a bit. We can chat more if you want, or skype $\endgroup$ Commented Oct 24, 2016 at 17:44
  • $\begingroup$ @tsuresuregusa many thanks for your comment , are you using it to compute capillary interactions between particles? Regarding the fem solution I know that the solution will depend on the mesh and on the length of the domain. However as you increase the size of the domain the solution converges to the one for unconfined particles as I checked with analytical solutions. Also the error in the integral should decrease increasing the numebrr of elements on the boindary (irrespectively of the problem you are simulating) but this is not the case. And this is the key point. $\endgroup$
    – SSC Napoli
    Commented Oct 24, 2016 at 18:29
  • $\begingroup$ I'm computing the diffusophoretic force between particles, so first need to solve the diffusion equation, get the tangential gradient of the concentration and use that as boundary condition on the stoke equations. You are correct saying that the integral should converge to the real value, the problem is that for getting accurate results you should decrease the elements in the whole domain, making it too slow to compute. I've worked several months trying to obtain accurate results but finally decided to change the method to one only dependent on the boundary. $\endgroup$ Commented Oct 24, 2016 at 19:45
  • $\begingroup$ from what I could see, the problem seems to be the what you plot in the last plot in yellow. The solution of fem is accurate only on the nodes, and it fits a 2nd order function in between, that when you derive is highly oscillatory. One thing I tried to do was to undersample the result in say 20-50 points and interpolate a spline through it, that at least gave a more continuous derivative. $\endgroup$ Commented Oct 24, 2016 at 19:48

1 Answer 1


More an extended comment than an answer but may help to illustrate my last comment.

If you plot the solution normalised by the boundary condition, you can see that making the mesh finer creates a high frequency noise. If you differentiate this, say with respect to the angle, will have very large values on the nodes of the elements. This is for nptspart=2000

Plot[ufun1[dist/2 + Cos[θ], 
   Sin[θ]]/(1 + (Cos[θ])^2 - 
    Sin[θ]^2), {θ, -.01, .02 π}]

enter image description here

I don't know the internals of how the derivative of an interpolated function is computed but I guess this noise is what makes the integral so time consuming.

As a workaround I found that interpolating in fewer points gives you a smoother solution. This is the code I use, which is not directly applicable to your problem but will give you an idea:

First I compute the gradient

gradC[r_, z_] := Evaluate[D[concentration[r, z], {{r, z}}]];

then the tangential component of the gradient on each particle (mine are in +- z0)

vSlip[r_, z_] := 
  Module[{ϕ = ArcTan[r, Abs[z] - z0]}, 
     z] - {r, z - Sign[z] z0} ({r, z - Sign[z] z0}.gradC[r, z])];

And with that create an interpolated function fixing the values so it goes to zero on the axis (this was were my problem was)

fInterpol1 = 
     If[Abs[ϕ] >= π/2, 0, 
      vSlip[Cos[ϕ], z0 + Sin[ϕ]][[
       1]]]}, {ϕ, -π, π, π/5}], Method -> "Spline"];
fInterpol2 = 
     If[Abs[ϕ] >= π/2, 0, 
      vSlip[Cos[ϕ], z0 + Sin[ϕ]][[
       2]]]}, {ϕ, -π, π, π/5}], Method -> "Spline"];

with this, I can use a more smooth function that's fast to integrate.

vSlip2[r_, z_] := 
  Module[{ϕ = ArcTan[r, Abs[z] - z0]}, {fInterpol1[ϕ], 

In your case you should interpolate in 2D as I think yours is the full gradient that is used.

However, we were interested on how the integral of vslip changed with the distance between the particles and due to the errors we could not find anything conclusive, thus we are using BEM now.

  • $\begingroup$ Many thanks again for posting this comment, i believe that this could potentially fix the integration. However, my concerns on the basics still remain; What is the cause of that high frequency noise? Is it possible to get rid of it? In principle on a boundary element the value of the solution should be smooth as it is interpolated quadratically between the closest mesh node values...My idea is that the interpolation function is being evaluated outside the domain where it is defined, but i don't know why i don't get any warning about that... $\endgroup$
    – SSC Napoli
    Commented Oct 26, 2016 at 9:05

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.