7
$\begingroup$

I am considering the following integral:

$$J(x_1^2, x_2^2) := \frac{3 \cdot 2^{15} \pi^4}{96} x_1^2 x_2^2 \int_{-\infty}^\infty d\tau_3 \int_{-\infty}^\infty d\tau_4 \int_{-\infty}^\infty d\tau_5 \int_{-\infty}^\infty d\tau_6\ \Theta(\tau_{34}) \Theta(\tau_{45}) \Theta(\tau_{56}) (I_{13} I_{25} + I_{15} I_{23})(I_{14} I_{26} + I_{16} I_{24}), \tag{1}$$

with $\tau_{ij} := \tau_i - \tau_j$, $\Theta(x)$ the Heaviside step function and with:

$$I_{ij} := \frac{1}{(2\pi)^2} \frac{1}{x_i^2 + \tau_j^2}. \tag{2}$$

The awkward numerical prefactor in $(1)$ is intended to normalize $J(x_1^2,x_2^2)$ such that:

$$\lim\limits_{x_2 \to 0} J(1, x_2^2) = 1. \tag{3}$$

Now my goal is the following: compute $J(1,x_2^2)$ for $x_2$ as close to $0$ as possible, and as accurately as possible. I can do the $\tau_4$ and $\tau_6$ integrals analytically, thus leaving a $2$-dimensional integral to be done numerically. The issue I have is the following: the integral tends towards $1$ until some value of $x_2$, then suddenly drops close to $0$ when I further decrease $x_2$. I can push to smaller values of $x_2$ by increasing the precision goal and MinRecursion, but even doing so I do not manage to go below $x_2 = 10^{-6}$, for which I obtain:

$$J(1, 10^{-6}) = (0.999988142 \pm .000000014). \tag{4}$$

The context: I am trying to solve numerically the following equation:

$$J(1, x_2 \to 0) = 1 + c \cdot x_2 + \mathcal{O}(x_2^2). \tag{5}$$

In order to obtain an accurate value for $c$, I need to reduce as much as possible the contamination from the quadratic term. I would appreciate reaching an accuracy of $10^{-10}$ or even better if possible.

My code so far:

I13 = 1/(2 π)^2 1/(x1^2 + τ3^2);
I14 = 1/(2 π)^2 1/(x1^2 + τ4^2);
I15 = 1/(2 π)^2 1/(x1^2 + τ5^2);
I16 = 1/(2 π)^2 1/(x1^2 + τ6^2);
I23 = 1/(2 π)^2 1/(x2^2 + τ3^2);
I24 = 1/(2 π)^2 1/(x2^2 + τ4^2);
I25 = 1/(2 π)^2 1/(x2^2 + τ5^2);
I26 = 1/(2 π)^2 1/(x2^2 + τ6^2);

Assuming[{x1 > 0, x2 > 0, τ3 ∈ Reals, τ4 ∈ Reals, τ5 ∈ Reals},
  Integrate[
    3*2^15 π^4 1/96 x1^2 x2^2 (I13*I25 + I15*I23)*(I14*I26 + I16*I24)*
     HeavisideTheta[τ5 - τ6] HeavisideTheta[τ4 - τ5] \
HeavisideTheta[τ3 - τ4] HeavisideTheta[τ3 - τ5], {τ6, -∞, ∞}]] // FullSimplify;
integrand2 = 
 Assuming[{x1 > 0, x2 > 0, τ3 ∈ Reals, τ4 ∈ Reals}, 
    Integrate[%, {τ4, -∞, ∞}]] // Normal // FullSimplify

x2 = 10^-6;
AG = 15;
PG = 15;
WP = 50;
NIntegrate[SetPrecision[integrand2 /. {x1 -> 1}, 
  WP + 2], {τ3, -∞, ∞}, {τ5, -∞, ∞}, 
 Method -> {"GlobalAdaptive", "SymbolicProcessing" -> 0, 
   "SingularityHandler" -> None, MinRecursion -> 9}, PrecisionGoal -> AG, 
 AccuracyGoal -> PG, WorkingPrecision -> WP, 
 IntegrationMonitor :> ((errors = Through[#1@"Error"]) &)]
Total@errors
Clear[x2]

Edit: following the advice in the comments, here is the result of the Integrate to save some time:

integrand2 = (1/(2048 π^8)) x1 x2 (1/((x2^2 + τ3^2) (x1^2 + τ5^2)) + 
    1/((x1^2 + τ3^2) (x2^2 + τ5^2))) (-π ArcTan[τ5/x1] + 
    ArcTan[τ3/x2] (π + 2 ArcTan[τ5/x1]) - (π + 4 ArcTan[τ5/x1]) ArcTan[τ5/x2] + 
    ArcTan[τ3/x1] (π + 2 ArcTan[τ5/x2])) HeavisideTheta[τ3 - τ5]

Edit 2: Following a discussion with yarchik, it should be made clear whether the limit $x_2 \to 0$ and the integrations can be commuted. Using a simpler example, we can see that the limit $x_2 \to 0$ causes problem, while sending $x_2$ to a finite non-zero value remains okay. Here is a script for the simpler example:

I12 = 1/(2 π)^2 1/(x1-x2)^2;
I13 = 1/(2 π)^2 1/(x1^2 + τ3^2;
I24 = 1/(2 π)^2 1/(x2^2 + τ4^2);
integrand=x2^2*I12*I13*I24;
D[integrand, x2];
Assuming[x1 > 0,Integrate[%, {τ3, -∞, ∞}, {τ4, -∞, ∞}]];
Limit[%, x2 -> 0, Direction -> "FromAbove"]

In that case, if you switch the limit and the integral you get an erroneous zero.

Improve this question
$\endgroup$
4
  • $\begingroup$ By the way, note that if I change $x_2$ to $10^{-5}$ in the code above, the error is reduced to $10^{-9}$, which makes me believe that the accuracy given in eq. $(4)$ is somewhat doubtful. $\endgroup$
    – Pxx
    Jun 7, 2020 at 12:02
  • $\begingroup$ A small update: with GaussKronrod I seem to be able to go below $10^{-6}$, and the error does follow the accuracy goal that I specify. How much can I trust this method and the results given? (at $10^{-10}$ it starts decreasing again though) $\endgroup$
    – Pxx
    Jun 7, 2020 at 15:00
  • $\begingroup$ After having tried your code, I'd say that a major hurdle in getting any help is the fact that it takes FOREVER to get an answer to the Integrate, which you need to set up the integrand. Since your problem seems to be with the numerical integration, providing the pre-calculated results of the symbolic step would go a long way towards not wasting your bounty prize. Dump its InputForm into a pastebin or some such; do NOT provide a nb file though; those contain executable code, and most people don't open executable files from unknown sources. $\endgroup$
    – MarcoB
    Jun 13, 2020 at 0:50
  • $\begingroup$ @MarcoB Done! That's strange, on my computer (which is far from being a fast horse) it takes a few seconds, something like 15 seconds. $\endgroup$
    – Pxx
    Jun 13, 2020 at 7:45

2 Answers 2

Reset to default
4
+25
$\begingroup$

One can use the Leibnitz integral rule in order to demonstrate that the constant $c$ in the asymptotic expansion with respect to $x_2$ $$ J(1,x_2)\stackrel{x_2\rightarrow0}{\rightarrow}1+c x_2 +\mathcal{O}(x_2^2) $$ is zero. We have $$ c=\left.\frac{\mathrm{d}J(1,x_2)}{\mathrm{d}x_2}\right|_{x_2=0}=\lim_{x_2\rightarrow0}\frac{\mathrm{d}}{\mathrm{d}x_2}\Big[\int_{-\infty}^\infty\!\! d\tau_3 \int_{-\infty}^\infty \!\!d\tau_4 \int_{-\infty}^\infty \!\!d\tau_5 \int_{-\infty}^\infty\! \!d\tau_6\, F[x_1,x_2;\tau_3,\tau_4,\tau_5,\tau_6]\Big]\\ =\int_{-\infty}^\infty\! d\tau_3 \int_{-\infty}^\infty\! d\tau_4 \int_{-\infty}^\infty\! d\tau_5 \int_{-\infty}^\infty \!d\tau_6\ \lim_{x_2\rightarrow0}\frac{\mathrm{d}}{\mathrm{d}x_2}F[x_1,x_2;\tau_3,\tau_4,\tau_5,\tau_6]. $$ With the help of MA we get for the derivative of the integrand $F[x_1,x_2;\tau_3,\tau_4,\tau_5,\tau_6]$:

(D[3*2^15 π^4 1/96 x1^2 x2^2 (I13*I25 + I15*I23) (I14*I26 + I16*I24), x2]) /. {x1 -> 1, x2 -> 0}
(*0*)

Notice also that $x_1$ and $x_2$ enter the integrand quadratically.

Edit

A question about the legitimacy of the Leibnitz rule has been asked. Indeed, in the second example the direct computation and by the Leibnitz rule differ. In particular we have

I12=1/(2 π)^2 1/(x1-x2)^2;
I13=1/(2 π)^2 1/(x1^2+τ3^2);
I24=1/(2 π)^2 1/(x2^2+τ4^2);
int=x2^2*I12*I13*I24;
r1=Integrate[int,{τ3,-∞,∞},{τ4,-∞,∞},Assumptions->x1>0&&x2>0]

(*x2/(64 π^4 x1 (x1-x2)^2)*)

Series[r1,{x2,0,2}]

(*x2/(64 π^4 x1^3)+x2^2/(32 π^4 x1^4)+O[x2]^3*)

Exchanging the order of derivative/integration we get zero

(D[int,x2])/.{x2->0}
(*0*)

However, the structure of this integral is slightly different, meaning that the integrand int needs to be regularized.

I24b=1/(2 π)^2 1/((x2+a)^2+τ4^2);
regint=(x2+a)^2 I12 I13 I24b;
r2=(D[regint,x2])/.{x2->0};
Integrate[r2,{τ3,-∞,∞},{τ4,-∞,∞},Assumptions->x1>0&&x2>0&&a>0]/.a->0

(*1/(64 π^4 x1^3)*)

Now the prefactor of $x_2$ in both ways is the same: $$ \frac{1}{64 \pi ^4 x_1^3}. $$

We can do the same regularization to the original integral, however, the result is zero as in the original solution.

Improve this answer
$\endgroup$
10
  • $\begingroup$ I am not sure I understand why you think that the integral is zero based on the derivative? You are saying that eq. $(4)$ in the OP is not correct? Looking at my project as.a whole, it seems to be a consistent value. $\endgroup$
    – Pxx
    Jun 13, 2020 at 22:09
  • $\begingroup$ I see, that's a great idea! Thanks! Is it actually clear that we can commute the limit and the integrals? I mean, I do it in my original post but your answer got me thinking about it. $\endgroup$
    – Pxx
    Jun 14, 2020 at 10:00
  • 1
    $\begingroup$ I now actually think that we can't, but I might be wrong. Here is a simpler example. Consider integrand=x2^2*I12*I13*I24 integrated over $\tau_3$ and $\tau_4$ (MA can handle that easily). Applying your method, we can exchange the derivative and the integral without problem, but interchanging the limit and the integral gives an erroneous $0$. You can try it out with this code: D[integrand, x2]; Assuming[x1 > 0,Integrate[%, {\[Tau]3, -\[Infinity], \[Infinity]}, {\[Tau]4, -\[Infinity], \[Infinity]}]]; Limit[%, x2 -> 0, Direction -> "FromAbove"] $\endgroup$
    – Pxx
    Jun 14, 2020 at 10:30
  • $\begingroup$ Note that the problem arises only at $0$. If you put $x_2 \to 10^{-6}$, then the limits can be interchanged so the script of the original question is still valid I guess. $\endgroup$
    – Pxx
    Jun 14, 2020 at 10:39
  • $\begingroup$ Yes of course, I will also include the definition of $I_{12}$. $\endgroup$
    – Pxx
    Jun 14, 2020 at 10:46
3
$\begingroup$

Mathematica uses the floating-point model for number representation. This inherently gets Integrate and NIntegrate specially the fame of not generating much provable decimals at all.

Division is not the place for Mathematicas fame overall. Maple for example compares open their product to Mathematica on a focus among others for recursion. Wolfram Inc. did react and upgraded with some functions dedicated to recursion.

The named problem still remains in Mathematica:

s[i_] := s[i] = 2*s[i - 1] - 3*s[i - 1]^2;

s[0] = SetAccuracy[3/10, 20]

0.3000000000000000000

{s[1], s[2], s[10], s[25], s[40]}

{0.3300000000000000000, 0.333300000000000000, 0.33333333333333, \
0.3333, 0.*10^62}

Enhancement for Mathematica from Wolfram Incs is:

RecurrenceTable[{a[n + 1] == 2 a[n] - 3*a[n]^2, a[0] == 3/10}, a, {n, 
  1, 10}]

{33/100, 3333/10000, 33333333/100000000, \
3333333333333333/10000000000000000, \
33333333333333333333333333333333/100000000000000000000000000000000, \
3333333333333333333333333333333333333333333333333333333333333333/\
10000000000000000000000000000000000000000000000000000000000000000, 
 333333333333333333333333333333333333333333333333333333333333333333333\
33333333333333333333333333333333333333333333333333333333333/
  10000000000000000000000000000000000000000000000000000000000000000000\
0000000000000000000000000000000000000000000000000000000000000, 
 333333333333333333333333333333333333333333333333333333333333333333333\
3333333333333333333333333333333333333333333333333333333333333333333333\
3333333333333333333333333333333333333333333333333333333333333333333333\
33333333333333333333333333333333333333333333333/
  10000000000000000000000000000000000000000000000000000000000000000000\
0000000000000000000000000000000000000000000000000000000000000000000000\
0000000000000000000000000000000000000000000000000000000000000000000000\
0000000000000000000000000000000000000000000000000, 
 333333333333333333333333333333333333333333333333333333333333333333333\
3333333333333333333333333333333333333333333333333333333333333333333333\
3333333333333333333333333333333333333333333333333333333333333333333333\
3333333333333333333333333333333333333333333333333333333333333333333333\
3333333333333333333333333333333333333333333333333333333333333333333333\
3333333333333333333333333333333333333333333333333333333333333333333333\
3333333333333333333333333333333333333333333333333333333333333333333333\
33333333333333333333333/
  10000000000000000000000000000000000000000000000000000000000000000000\
0000000000000000000000000000000000000000000000000000000000000000000000\
0000000000000000000000000000000000000000000000000000000000000000000000\
0000000000000000000000000000000000000000000000000000000000000000000000\
0000000000000000000000000000000000000000000000000000000000000000000000\
0000000000000000000000000000000000000000000000000000000000000000000000\
0000000000000000000000000000000000000000000000000000000000000000000000\
0000000000000000000000000, 
 333333333333333333333333333333333333333333333333333333333333333333333\
3333333333333333333333333333333333333333333333333333333333333333333333\
3333333333333333333333333333333333333333333333333333333333333333333333\
3333333333333333333333333333333333333333333333333333333333333333333333\
3333333333333333333333333333333333333333333333333333333333333333333333\
3333333333333333333333333333333333333333333333333333333333333333333333\
3333333333333333333333333333333333333333333333333333333333333333333333\
3333333333333333333333333333333333333333333333333333333333333333333333\
3333333333333333333333333333333333333333333333333333333333333333333333\
3333333333333333333333333333333333333333333333333333333333333333333333\
3333333333333333333333333333333333333333333333333333333333333333333333\
3333333333333333333333333333333333333333333333333333333333333333333333\
3333333333333333333333333333333333333333333333333333333333333333333333\
3333333333333333333333333333333333333333333333333333333333333333333333\
333333333333333333333333333333333333333333333/
  10000000000000000000000000000000000000000000000000000000000000000000\
0000000000000000000000000000000000000000000000000000000000000000000000\
0000000000000000000000000000000000000000000000000000000000000000000000\
0000000000000000000000000000000000000000000000000000000000000000000000\
0000000000000000000000000000000000000000000000000000000000000000000000\
0000000000000000000000000000000000000000000000000000000000000000000000\
0000000000000000000000000000000000000000000000000000000000000000000000\
0000000000000000000000000000000000000000000000000000000000000000000000\
0000000000000000000000000000000000000000000000000000000000000000000000\
0000000000000000000000000000000000000000000000000000000000000000000000\
0000000000000000000000000000000000000000000000000000000000000000000000\
0000000000000000000000000000000000000000000000000000000000000000000000\
0000000000000000000000000000000000000000000000000000000000000000000000\
0000000000000000000000000000000000000000000000000000000000000000000000\
00000000000000000000000000000000000000000000000}

This rather short excursion sheets some light on the kernel of the given question.

New in version 12 is the built-in Typed. It introduces a representation of expressions that should be assumed to be of a specified type for compilation and other purposes.

An example from FunctionCompile is:

cfun = FunctionCompile[
  Function[{Typed[arg, "MachineInteger"]}, 
   Module[{factFun}, 
    factFun = 
     Function[{arg1}, If[arg1 === 1, 1, arg1*factFun[arg1 - 1]]];
    factFun[arg]]]];

First[RepeatedTiming[cfun[10]]/RepeatedTiming[10!]]

5.86162269785962781826373415267759306657260113508390167035153045127714\
38976022830

This throws warning concerning General::precsm and gives more precision comparable to Maple or Matlab or C/C++.

So for more precision put all the functions into the FunctionCompile built-in and invest the time.

I agree with the suggestion from the other answer so far it advances the precision of the result.

Some enhancement can be taken from failed to use setprecision

The example is:

ER2[α_, K_, q_] := 
  1 + Sum[Binomial[α + K, m] Sum[
      Binomial[m, r]*(-1)^r/((1/q)^(α + K + r - m) - 1), {r, 0,
        m}], {m, 0, K - 1}];
ER2[0, 268, 1/10] // N
ER2[0, 268, SetPrecision[0.1, 50]]
SetPrecision[Unevaluated@ER2[0, 268, 0.1], 50]


3.19697
3.196969869247707935345524749
3.1969698692477079353455247

But that is built-in Mathematica precision enhancement.

Another chance is the use of NSum.

{sum1, {points1}} = 
  Reap[NSum[Log[i^2]/(2^i i!), {i, ∞}, 
    Method -> "WynnEpsilon", EvaluationMonitor :> Sow[i]]];
{sum2, {points2}} = 
  Reap[NSum[Log[i^2]/(2^i i!), {i, ∞}, 
    Method -> "EulerMaclaurin", EvaluationMonitor :> Sow[i]]];
{sum1 - sum2, points1 // Length, points2 // Length}

{ListPlot[Transpose@{points1, Range[Length[points1]]}, 
  ImageSize -> 200], 
 ListPlot[Transpose@{points2, Range[Length[points2]]}, 
  ImageSize -> 200]}

ListPlot of convergence

ListPlot of convergence

This is an example out of the book Mathematica Navigator, page 669. It gives more insight into the convergence of the series and adds some theory about solving methodology for such kind of series.

Shows what is written there about the methods.

Both methods offer precision control with PrecisionGoal and AccuracyGoal. The Mathematica documentation lists the WynnEpsilon under the section "Possible Issues".

Wynn's extrapolation algorithm loses precision for series that are not alternating:

NSum[1/n^3, {n, 1, ∞}, Method -> "WynnEpsilon", 
  WorkingPrecision -> 30] // Precision

7.27779

The series in this question appear not to be alternating.

The The Mathematica GuideBook for Numericsgives an exmaple for this question on page 281:

Timing[NSum[1/(i^2 + 2), {i, Infinity}, Method -> #]] & /@ {Integrate,
    NIntegrate, SequenceLimit, Fit} // InputForm

{{0.296279, 0.8610281004917035}, {0.004104, 0.8610281004917042}, 
 {0.001149, 0.8598255911884887}, {0.001189, 0.8598255911884887}}

The terms are similar to the ones of this question. Mr. Michael Trott recommends therefore the default method as already chosen and find an optimum with the EvaluationMonitor method already presented.

Mr. Trott presents in that book examples in which VerifyConvergence->False and Off[SequenceLimit::seqlim] work. Have a try for advances in more precision. All can be ready in the Google books online version of the book freely available. I was not able to find SequenceLimit in the Mathematica documentation but there is a built-in named NLimit with the option SequenceLimit. But the examples work in 12.0 for me.

That is it from the Mathematica inside for this type of numerical calculations. Always keep in mind that from the university standpoint a remainder calculation is necessary that leads to the O-notation calculus.

The preliminaries thoughs are needed to the active for reading the article:

NIntegrate Integration Strategies.

Since the strategies are:

"The "LocalAdaptive" strategy has an initialization routine and a Recursive Routine (RR). RR produces the leaves of a tree, the nodes of which are regions. The children of a node/region are subregions obtained by its partition. RR takes a region as an argument and returns an integral estimate for it."

as stated in the article from Wolfram Inc. in the section "Local Adaptive Strategy".

In this article, there is a section dealing with the integration of

Plot[1/(10 (1/2 - x)^2 + 1), {x, 0, 1}]

Plot

It is than demonstrated that:

"Multidimensional example of using the "Partitioning" option. To make the plot, the sampling points of the first region to be integrated, {0,1}x{0,1}, are removed:"

sampledPoints = 
  Reap[NIntegrate[(1 + (x + y)^2)^-1, {x, 0, 1}, {y, 0, 1}, 
     Method -> {"LocalAdaptive", "Partitioning" -> {3, 4}}, 
     EvaluationMonitor :> Sow[{x, y}]]][[2, 1]];
sampledPoints = 
  Partition[sampledPoints, 17(*Length[sampledPoints]/(3*4+1)*)];
sampledPoints = Flatten[Rest[sampledPoints], 1];
ListPlot[sampledPoints, AspectRatio -> 1, 
 GridLines -> {Range[3]/3, Range[4]/4}]

Sampeled points

It follows that for functions of the type 1/(a+(x+y)^2) or alike are well integrated with the "Example Implementation of a Local Adaptive Strategy" from the referred Wolfram Inc. article. Nice precision goals are easily and fast and with high reproducibility achieved.

This table extends the answer:

"GlobalAdaptive" vs. "LocalAdaptive"

I recommend for the numerical part to evaluate the article in full perhaps some nice ideas to get more precision will be gained in the special case, but the article focusses on singularities more than 2-dimensional integration precision. It shows up many test on convergence and precision that are fun to apply. The article has 85 pages in length.

Here are calculation for reference with the Mathematica example for WorkingPrecision:

AG = 15;
PG = 15;
WP = 50;
Table[NIntegrate[SetPrecision[1/Sqrt[Sin[x]], wp + 2], {x, 0, 1}, 
   Method -> {"GlobalAdaptive", "SymbolicProcessing" -> 0, 
     "SingularityHandler" -> None, MinRecursion -> 9}, 
   PrecisionGoal -> AG, AccuracyGoal -> PG, WorkingPrecision -> WP, 
   IntegrationMonitor :> ((errors = Through[#1@"Error"]) &)] - 
  integexact, {wp, 20, 40, 2}]
Total@errors

{-0.0000403805375729249773658836249940223892505415071 + (
  1176456130984411 I)/
  10141204801825835211973625643008, \
-0.0000403805375729249773658836249940223892505415071 + (
  1176456130984411 I)/
  10141204801825835211973625643008, \
-0.0000403805375729249773658836249940223892505415071 + (
  1176456130984411 I)/
  10141204801825835211973625643008, \
-0.0000403805375729249773658836249940223892505415071 + (
  1176456130984411 I)/
  10141204801825835211973625643008, \
-0.0000403805375729249773658836249940223892505415071 + (
  1176456130984411 I)/
  10141204801825835211973625643008, \
-0.0000403805375729249773658836249940223892505415071 + (
  1176456130984411 I)/
  10141204801825835211973625643008, \
-0.0000403805375729249773658836249940223892505415071 + (
  1176456130984411 I)/
  10141204801825835211973625643008, \
-0.0000403805375729249773658836249940223892505415071 + (
  1176456130984411 I)/
  10141204801825835211973625643008, \
-0.0000403805375729249773658836249940223892505415071 + (
  1176456130984411 I)/
  10141204801825835211973625643008, \
-0.0000403805375729249773658836249940223892505415071 + (
  1176456130984411 I)/
  10141204801825835211973625643008, \
-0.0000403805375729249773658836249940223892505415071 + (
  1176456130984411 I)/10141204801825835211973625643008}


0.000062690729178774772783412492441137149675804414387926

output of ListPlot for the data

Table[NIntegrate[SetPrecision[1/Sqrt[Sin[x]], wp + 2], {x, 0, 1}, 
   PrecisionGoal -> AG, AccuracyGoal -> PG, WorkingPrecision -> wp, 
   IntegrationMonitor :> ((errors = Through[#1@"Error"]) &)] - 
  integexact, {wp, 20, 40, 2}]
Total@errors

{4.972*10^-16 + (1176456130984411 I)/10141204801825835211973625643008,
  4.97155*10^-16 + (1176456130984411 I)/
  10141204801825835211973625643008, 
 4.9715456*10^-16 + (1176456130984411 I)/
  10141204801825835211973625643008, 
 4.971545601*10^-16 + (1176456130984411 I)/
  10141204801825835211973625643008, 
 4.97154560077*10^-16 + (1176456130984411 I)/
  10141204801825835211973625643008, 
 4.9715456007147*10^-16 + (1176456130984411 I)/
  10141204801825835211973625643008, 
 4.971545600714668*10^-16 + (1176456130984411 I)/
  10141204801825835211973625643008, 
 4.97154560070974729*10^-16 + (1176456130984411 I)/
  10141204801825835211973625643008, 
 4.9715456007097409096*10^-16 + (1176456130984411 I)/
  10141204801825835211973625643008, 
 4.971545600709740909560*10^-16 + (1176456130984411 I)/
  10141204801825835211973625643008, 
 4.97154560070974083528010*10^-16 + (1176456130984411 I)/
  10141204801825835211973625643008}

4.93444152315623988446577100475728722424000980*10^-16

Happy use of Mathematica!

Improve this answer
$\endgroup$
5
  • $\begingroup$ That is very interesting, as even with the solution suggested by yarchik I will have to do other similar computations which will not always be zero. I will try it out and report what comes out. $\endgroup$
    – Pxx
    Jun 14, 2020 at 21:42
  • $\begingroup$ If I copy you code in Mathematica, I always get a "beep" and results such as '9.', '4.' or '3.', any idea why? I am using version 12.0.0.0. $\endgroup$
    – Pxx
    Jun 16, 2020 at 11:44
  • $\begingroup$ I addressed the example you open with in (60124) and I still do not see it as a weakness of Mathematica. $\endgroup$
    – Mr.Wizard
    Jun 17, 2020 at 8:54
  • $\begingroup$ @Mr.Wizard I actually was able to increase the precision of the results by setting a higher precision of the variable $x_2$ (in the same way as you suggest in the post that you mention), but the results appear to be non-sensical: the coefficient $c$ seems to diverges as I approach $0$! $\endgroup$
    – Pxx
    Jun 17, 2020 at 9:08
  • $\begingroup$ Unfortunately, I cannot give the bounty to this answer. I have tried everything here (except the NSum as I did not know how to do it), and most of the suggested solutions (e.g. using LocalAdaptive) led to the same problem or did not work. So far the only thing that increased the range for my integral was to use SetPrecision, but the results are still not satisfying. Thanks for your contribution though, still I have learnt a great deal! $\endgroup$
    – Pxx
    Jun 18, 2020 at 7:54

Your Answer

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

Not the answer you're looking for? Browse other questions tagged or ask your own question.