7
$\begingroup$

I have mathematica 14.3.

I have encountered a weird issue when numerically integrating this complicated function

f[x_] := 
  UnitStep[
     x - .3]*(1/(-0.3` + x)^2 (-1.1434102576408818` - 
        73.9413929141465` (-0.3` + x)^2.915` - 
        115.82390213796239` (-0.3` + x)^3.415` + 
        536.3471523421645` (-0.3` + x)^3.915` - 
        513.5287667970828` (-0.3` + x)^4.415` + 
        24.579421179117183` (0.3` + x)^2.915` + 
        24.925561694460875` (0.3` + x)^3.415` - 
        77.76754995587882` (0.3` + x)^3.915` + 
        51.33230871908811` (0.3` + x)^4.415` + 
        x (2.65606615372176` - 
           154.94956890020973` (-0.3` + x)^2.915` - 
           145.89473463682685` (-0.3` + x)^3.415` + 
           388.39637357704294` (-0.3` + x)^3.915` - 
           178.45013745326605` (-0.3` + x)^4.415` - 
           124.82340431573319` (0.3` + x)^2.915` - 
           114.87401830151424` (0.3` + x)^3.415` + 
           331.86338927447434` (0.3` + x)^3.915` - 
           205.70188178134828` (0.3` + x)^4.415`) + 
        x^2 (-167.90881183998138` (-0.3` + x)^2.915` - 
           145.10636075108883` (-0.3` + x)^3.415` + 
           412.2021927513278` (-0.3` + x)^3.915` - 
           266.5287241054068` (-0.3` + x)^4.415` + 
           165.54009634792362` (0.3` + x)^2.915` + 
           139.07326720934407` (0.3` + x)^3.415` - 
           372.5592895364412` (0.3` + x)^3.915` + 
           216.691879843441` (0.3` + x)^4.415`))) + 
   UnitStep[.3 - 
      x]*(1/(-0.3` + x)^2 (-1.1434102576408818` + 
        24.579421179117187` (0.3` - x)^2.915` + 
        24.925561694460878` (0.3` - x)^3.415` - 
        77.76754995587882` (0.3` - x)^3.915` + 
        51.33230871908812` (0.3` - x)^4.415` + 
        24.579421179117187` (0.3` + x)^2.915` + 
        24.925561694460878` (0.3` + x)^3.415` - 
        77.76754995587882` (0.3` + x)^3.915` + 
        51.33230871908812` (0.3` + x)^4.415` + 
        x (2.65606615372176` + 
           124.82340431573319` (0.3` - x)^2.915` + 
           114.87401830151425` (0.3` - x)^3.415` - 
           331.86338927447434` (0.3` - x)^3.915` + 
           205.70188178134828` (0.3` - x)^4.415` - 
           124.82340431573319` (0.3` + x)^2.915` - 
           114.87401830151425` (0.3` + x)^3.415` + 
           331.86338927447434` (0.3` + x)^3.915` - 
           205.70188178134828` (0.3` + x)^4.415`) + 
        x^2 (165.54009634792365` (0.3` - x)^2.915` + 
           139.07326720934404` (0.3` - x)^3.415` - 
           372.5592895364412` (0.3` - x)^3.915` + 
           216.691879843441` (0.3` - x)^4.415` + 
           165.54009634792365` (0.3` + x)^2.915` + 
           139.07326720934404` (0.3` + x)^3.415` - 
           372.5592895364412` (0.3` + x)^3.915` + 
           216.691879843441` (0.3` + x)^4.415`)));

The weird thing is that

NIntegrate[f[x], {x, 0, 1}]

gives

enter image description here

-1162.46 - 5.61143*10^-32 I

while

NIntegrate[f[x], {x, 0, .4}] + NIntegrate[f[x], {x, .4, 1}]

gives the correct answer

7.5341

Why does this happen? Is there an easy fix without splitting up the integral at some random point by trial and error?

Improve this question
$\endgroup$
9
  • 2
    $\begingroup$ NIntegrate[f[x], {x, 0, 1}, MinRecursion -> 1] seems to produce an accurate result $\endgroup$ Commented Sep 15 at 22:17
  • 1
    $\begingroup$ If you reduce the accuracy goal it also works. NIntegrate[f[x],{x,0,1},AccuracyGoal->4] gives 7.5341. Btw, the point .3 is problem in your input. Limit[Rationalize[f[x]],x->3/10] gives Indeterminate same as for f[.3]. This the point where the unit steps join. But Exclusions -> {.3} did not help. Do not know if this is why you get the original error. i.sstatic.net/Xx2PQMcg.png $\endgroup$ Commented Sep 15 at 22:40
  • 2
    $\begingroup$ Another workaround: fe = FunctionInterpolation[f[x], {x, 0, 1}]; NIntegrate[fe[x], {x, 0, 1}] gives 7.53412 $\endgroup$ Commented Sep 15 at 23:01
  • 1
    $\begingroup$ The error explains the problem doesn't it? $\endgroup$ Commented Sep 15 at 23:32
  • 2
    $\begingroup$ @Nasser The error points to the numerical problem I show in my answer. In this case, if the numerical badness of the singularity is missed by NIntegrate[] (as happens in the integral over {x, 0, 0.4}), the error introduced by the "mistake" is negligible, because the "exact" integrand is actually well-behaved at the singularity. But if the badness is caught, then NIntegate[] crumbles under the numerical instability. The magnification of round-off error by the factors 1/(-0.3 + x)^2 for x near 0.3 ought to be considered in attempting to integrate f[x]. $\endgroup$ Commented Sep 16 at 1:21

2 Answers 2

Reset to default
6
$\begingroup$

It's numerics:

Plot[f[x], {x, 0.2999999999999, 0.3000000000001}]

plot with a range exceeding 10^12

The pink highlighting is because the plot generated a system exception error. Trying to Rasterize[] the output generates SystemException["MemoryAllocationFailure"]. I cannot tell what system exception the plot generated.

Note how large the values of the integrand are.

Now if NIntegrate[] succeeds, it's because the error estimator did not sample the unstable region. Or at least I assume so, because if it samples a few points in the above region, the estimated error will be very large and NIntegrate[] will eventually complain. It's not mysterious. If it misses the bad region, then the error estimate might be low enough for it to report success.

In any case, there are a couple of ways that seem to work besides dividing the region at magic places. One is to use high precision to tame the numerical instability. The other is to turn off symbolic processing, but I think that method might be akin to the magic-place approach: It happens to work this time.

NIntegrate[f[x] // Rationalize // Rationalize[#, 0] &, {x, 0, 3/10, 1},
  WorkingPrecision -> 60, PrecisionGoal -> 6]
(*  7.53410157432435900859372569223172462580179350407091264752486  *)

(* PrecisionGoal -> 8 fails if WorkingPrecision -> 60 *)
NIntegrate[f[x] // Rationalize // Rationalize[#, 0] &, {x, 0, 3/10, 1},
  WorkingPrecision -> 80, PrecisionGoal -> 8]
(*  7.5341016185204640407847706393718283467283168485844385359272221004580751170518902  *)

NIntegrate[f[x] // Rationalize // Rationalize[#, 0] &, {x, 0, 1}, 
 Method -> {"GaussKronrodRule", "SymbolicProcessing" -> None}]
(*  7.5341  *)

The failure of PrecisionGoal -> 8 at WorkingPrecision -> 60 suggests the numerical problems are significant.

Improve this answer
$\endgroup$
2
  • $\begingroup$ This shows the numerical noise at different scales around 0.3: Table[ListLinePlot[f[Subdivide[0.3 + 0.2*10^-k, 0.3 + 10^-k, 10000]] // RealExponent // {MaxFilter[#, 20], MinFilter[#, 20]} &, PlotRange -> All, PlotStyle -> Thickness@Small, DataRange -> {0.2*10^-k, 10^-k}], {k, 5, 7}] // GraphicsRow. Out around $0.3+10^{-4}$, the noise becomes less than $10^{-7}$, probably still large enough to bother NIntegrate[] but too small to see with the code. $\endgroup$ Commented Sep 16 at 11:07
  • $\begingroup$ Replacing UnitStep -> HeavisideTheta solves the problems. See my solution. $\endgroup$ Commented Sep 19 at 20:39
2
$\begingroup$

Define a function g[x] exactly like f[x] but with UnitStep replaced by HeavisideTheta and everything works out fine for g[x].

No problems with splitting the integration interval or accuracy.

f and g look identical in the plot

Plot[{g[x], f[x]}, {x, -1, 2}]

enter image description here

But the integral over f and g differs appreciately

In[97]:= NIntegrate[{g[x], f[x]}, {x, 0, 1}]

Out[97]= {7.5341, -1162.46 - 5.61143*10^-32 I}

and only that over g is correct.

Improve this answer
$\endgroup$

Your Answer

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

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.