2
$\begingroup$

I'm trying to integrate an interpolating function, but when I use NIntegrate, Mathematica gives me a series of error messages. I give you some examples:

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::inumri: The integrand <<1>> has evaluated to Overflow, Indeterminate, or Infinity for all sampling points in the region with boundaries {{14.6836,14.6953},{-12.0000170249736609318194678966418287302531098248437047004699707031,-12.0133777475052196583721331180072411370929330587387084960937500000}}.

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.

General::stop: Further output of NIntegrate::inumri will be suppressed during this calculation.

I hope someone will help me to solve this problem. Here is the code of a simplified version of my problem:

Ω = 
ImplicitRegion[
0 <= x <= 34 && -18 <= y <= 
0 && ! (12 < x < 15 && -12 < y || 15 < x < 20 && -3 < y || 
20 < x < 23 && -7 < y || 23 < x <= 34 && -3 < y), {x, y}];
RegionPlot[Ω, AspectRatio -> 0.5]

fi = NDSolve[{\!\(
\*SubsuperscriptBox[\(∇\), \({x, y}\), \(2\)]\(φ[x,
y]\)\) == NeumannValue[0, 
x == 0 && -18 < y < 0 || 0 < x < 34 && y == -18 || 
x == 34 && -18 < y < -3 || x == 23 && -7 < y < -3 || 
20 < x < 23 && y == -7 || x == 20 && -7 < y < -3 || 
15 < x < 20 && y == -3 || x == 15 && -12 < y < -3 || 
12 < x < 15 && y == -12 || x == 12 && -12 < y < 0], 
DirichletCondition[φ[x, y] == 1.5, 
0 <= x <= 12 && y == 0], 
DirichletCondition[φ[x, y] == 0.4, 
23 <= x <= 34 && y == -3]}, φ, {x, 
y} ∈ Ω]
φ[x_, y_] = φ[x, y] /. fi;

Integrale = 
NIntegrate[D[φ[x, y], x], {x, y} ∈ Ω]
$\endgroup$
1
  • $\begingroup$ This \[CapitalOmega] =ImplicitRegion[..YourStuff..]; fi=\[CurlyPhi][x,y] /. NDSolve[..YourStuff..][[1,1]]; Print[Plot3D[fi, {x,0,34}, {y,-18,0}]]; NIntegrate[D[fi,x], {x,y} \[Element] \[CapitalOmega]] fixes part of your problems with NIntegrate immediately failing, but still doesn't fix your convergence problem. $\endgroup$
    – Bill
    Commented Jan 6, 2017 at 1:32

2 Answers 2

5
$\begingroup$

First of all, on an InterpolatingFunction constructed by NDSolve one should use its ElementMesh[..] domain in NIntegrate, because that will invoke a special FEM integrator that is fast and accurate on the interpolating function.

One finds sometimes that the discretization of a region by NDSolve and by NIntegrate are not exactly the same, which can lead to sampling outside the domain. The problem can usually be fixed by using the ElementMesh created by NDSolve, but for some reason that did not work. (The integral evaluates to 0., which does not agree with the graph of the integrand, which is negative.) Further, sampling outside the domain seems still to occur and cause overflow and for NIntegrate to return unevaluated.

The next natural thing to try was to extend the function by zero outside the domain. One way this can be done with the option "ExtrapolationHandler". This works.

Experimenting with a coding way to extend the function by zero if the input is not in the region tended to be unsuccessful.

Clear[φ];
fi = NDSolve[{
   Derivative[0, 2][φ][x, y] + Derivative[2, 0][φ][x, y] == NeumannValue[0, 
     x == 0 && -18 < y < 0 || 0 < x < 34 && y == -18 || 
      x == 34 && -18 < y < -3 || x == 23 && -7 < y < -3 || 
      20 < x < 23 && y == -7 || x == 20 && -7 < y < -3 || 
      15 < x < 20 && y == -3 || x == 15 && -12 < y < -3 || 
      12 < x < 15 && y == -12 || x == 12 && -12 < y < 0], 
   DirichletCondition[φ[x, y] == 1.5, 0 <= x <= 12 && y == 0], 
   DirichletCondition[φ[x, y] == 0.4, 23 <= x <= 34 && y == -3]},
   φ, {x, y} ∈ Ω,
  "ExtrapolationHandler" -> {0 &, "WarningMessage" -> False}];

emesh = φ["ElementMesh"] /. First[fi];
NIntegrate[D[φ[x, y], x] /. First[fi], {x, y} ∈ emesh]
(*  -6.00409  *)

Alternatives

One needs to reevaluate the OP's NDSolve without the "ExtrapolationHandler" to test the following.

ClearAll[ff];
dxφ = Derivative[1, 0][φ] /. First[fi];
ff[x_?NumericQ, y_?NumericQ] := Check[dxφ[x, y], 0.];
NIntegrate[ff[x, y], {x, y} ∈ emesh]
(*  -6.00409  *)

These, all of which use emesh, return zero (all return unevaluated as in the OP if emesh is replaced by Ω):

NIntegrate[D[φ[x, y], x] /. First[fi], {x, y} ∈ emesh]
NIntegrate[dxφ[x, y], {x, y} ∈ emesh]
NIntegrate[Check[dxφ[x, y], 0.], {x, y} ∈ emesh]

The following works but with limitations. NIntegrate will tell you that you cannot use "GlobalAdaptive" on regions, but the error message suggests that it is in fact used as a submethod. But how to pass the submethod options? That's probably worth it's own Q&A.

ClearAll[ff];
ff[x_?NumericQ, y_?NumericQ] := Quiet@Check[dxφ[x, y], 0.];
NIntegrate[ff[x, y], {x, y} ∈ Ω]

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 2000 times....

    (*  -6.00189  *)
$\endgroup$
1
  • $\begingroup$ Very informative! $\endgroup$ Commented Jan 8, 2017 at 16:17
0
$\begingroup$

It works by defining the ImplicitRegion in a different way:

\[CapitalOmega]1 = ImplicitRegion[0 <= x <= 12 && -18 <= y <= 0, {x, y}];
\[CapitalOmega]2 = ImplicitRegion[12 <= x <= 15 && -18 <= y <= -12, {x, y}];
\[CapitalOmega]3 = ImplicitRegion[15 <= x <= 20 && -18 <= y <= -3, {x, y}];
\[CapitalOmega]4 = ImplicitRegion[20 <= x <= 23 && -18 <= y <= -7, {x, y}];
\[CapitalOmega]5 = ImplicitRegion[23 <= x <= 34 && -18 <= y <= -3, {x, y}];
\[CapitalOmega] = 
RegionUnion[\[CapitalOmega]1, \[CapitalOmega]2, \[CapitalOmega]3, \
\[CapitalOmega]4, \[CapitalOmega]5];

The rest of the code is the same.

$\endgroup$
1

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.