# Numerical integration converging too slowly

I must solve this integral which I suppose to be a very small number. How can I do? When I wrote this code:

r = (1.082)*10^8
h = (4.87)*10^(24)
Q[x_] := x - (0.206)*Sin[x] + (0.206)*(0.206)*Sin[x]*Cos[x]
R[x_] := ((57.91)*10^6)*(1 + 0.206*Cos[Q[x]])
B[x_] := (0.206)*Sin[2*x] - Sin[x]
A[x_] := ((0.62)*10^(-25))*((4.854)*Cos[Q[x]] - B[Q[x]]*Sin[Q[x]])
F[x_, y_] :=
R[x]*A[x] -
r*Cos[y]*((((57.91)*10^6)*A[x]*(1 - (11.93)*10^6)*
Cos[Q[x]]/R[x]) - ((-3.008)*10^(-25)) - ((57.91)*10^6)*
B[x]*((-0.052)*10^(-31))*Sin[Q[x]]) -
r*Sin[y]*(((0.6329)^10^(-25))*
Sin[Q[x]] + ((56.67)*10^6)*((-0.05196)*10^(-31))*B[x]*Cos[Q[x]])
G[x_, y_] := (1 - (2*
R[x]/r)*(Cos[y]*(((11.93)*10^6) + ((57.91)*10^6)*Cos[Q[x]])/
R[x]) + (R[x]/r)^2)^(-1.5)
u = ((0.01917)*10^(27))*((57.91)*10^6)
NIntegrate[(h/(((1.989)*10^(30)))*(r^3))*F[x, y]*G[x, y]*u, {x, 0,
2 Pi}, {y, 0, 2 Pi}, PrecisionGoal -> 50, MaxRecursion -> 50,
WorkingPrecision -> 100]


it gives me this error:**

"The precision of the argument function \
((3.44312*10^51\(<<1>>))/(1+0.286452\(<<1>>)^2-1.84843*10^-8\Cos[y]\(\
1.193*10^7+5.791*10^7\Cos[Plus[<<3>>]]))^1.5) is less than \
WorkingPrecision (100.)"


and this:

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

he global error of the strategy GlobalAdaptive has increased more \
than 2000 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 8.\
2245804667122205619329698365928375872918387045667313082675253184241492\
7464246649828330519930578973435118997853216575354898819456603870372689\
813756795*10^42 and 3.\
8524159983189941460321064133370663169006624537848246120142620883724145\
3794096694780053810840247925397951565645812529096240278267666686786349\
959791842*10^48 for the integral and error estimates.

• The error explicitly suggests using higher WorkingPrecision. Have you tried doing just that? In particular, my understanding is that you must set WorkingPrecision to be at least as large as PrecisionGoal. Set WorkingPrecision -> 1000? Commented Jun 11, 2015 at 16:54
• i set PrecisionGoal -> 50, MaxRecursion -> 50, WorkingPrecision -> 100 Commented Jun 11, 2015 at 17:10
• and he told me The precision of the argument function... is less than WorkingPrecision (100.) Commented Jun 11, 2015 at 17:11
• i tried also PrecisionGoal -> 1000, MaxRecursion -> 1000, WorkingPrecision -> 2000 Commented Jun 11, 2015 at 17:12
• OK, so what are the argument functions? Could you show how you defined them? Commented Jun 11, 2015 at 17:12

As noted in the Comments, this integration is plagued by precision problems. To proceed, factor the huge constant

(h/(((1.989)*10^(30)))*(r^3)) u
(* 3.44312*10^51 *)


from the integrand and then FullSimplify and Rationalize the functions F[x, y]*G[x, y]]

rat = Rationalize[FullSimplify[F[x, y]*G[x, y]], 0];


in order to achieve any possible cancellations symbolically (not many) and to increase the precision of the integrand to infinity. (The first error message in the Question was Mathematica complaining that the precision of the integrand was less (in fact, much less) than the requested WorkingPrecision -> 100.)

To gain some insight into the integrand,

Plot3D[rat, {x, 0, 2 Pi}, {y, 0, 2 Pi}, PlotRange -> All]


Again as noted in the comments, the integrand is highly symmetric in both dimensions. This suggests that the four quadrants be summed symbolically before integration.

Plot3D[(rat /. {x -> xx, y -> yy}) + (rat /. {x -> 2 Pi - xx,
y -> yy}) + (rat /. {x -> xx,
y -> 2 Pi - yy}) + (rat /. {x -> 2 Pi - xx,
y -> 2 Pi - yy}), {xx, 0,  Pi}, {yy, 0, Pi}, PlotRange -> All]


Observe that the modified integrand is some 15 orders of magnitude smaller (a huge cancellation) and also that it is highly irregular. It now can be integrated.

NIntegrate[(rat /. {x -> xx, y -> yy}) + (rat /. {x -> 2 Pi - xx,
y -> yy}) + (rat /. {x -> xx, y -> 2 Pi - yy}) + (rat /. {x -> 2 Pi - xx,
y -> 2 Pi - yy}), {xx, 0,  Pi}, {yy, 0, Pi},
MaxPoints -> 500000, WorkingPrecision -> 30, PrecisionGoal -> 9]
(* 3.94313*10^-9 *)


without error messages. Note that MaxPoints -> 500000, WorkingPrecision -> 30 is necessary to obtain reasonable accuracy, at the cost of a few minutes of computational time. Multiplying the answer by the huge constant factored out earlier gives 1.35767*10^43.

Improved computation

The question below prompted me to consider more carefully the second plot above and the subsequent integration. Applying FullSimplify to the argument of the second plot reduces its LeafCount from 887 to 192, and provides a much more accurate plot.

sim = FullSimplify[(rat /. {x -> xx, y -> yy}) + (rat /. {x -> 2 Pi - xx, y -> yy})
+ (rat /. {x -> xx, y -> 2 Pi - yy}) + (rat /. {x -> 2 Pi - xx, y -> 2 Pi - yy})];
Plot3D[sim, {xx, 0, Pi}, {yy, 0, Pi}, PlotRange -> All]


Then,

NIntegrate[sim, {xx, 0, Pi}, {yy, 0, Pi}, MaxPoints -> 1000000,
WorkingPrecision -> 30, PrecisionGoal -> 9]
(* -2.25620530195937290197090585866*10^-18 *)


Multiplying the answer by the huge constant factored out earlier gives -7.76838*10^33.

• Anyone kindly explain this " Plot3D[(rat /. {x -> xx, y -> yy}) + (rat /. {x -> 2 Pi - xx, y -> yy}) + (rat /. {x -> xx, y -> 2 Pi - yy}) + (rat /. {x -> 2 Pi - xx, y -> 2 Pi - yy}), {xx, 0, Pi}, {yy, 0, Pi}, PlotRange -> All] " .Please if possible kindly explain what is happening here? Commented Nov 25, 2019 at 13:01
• @GummalaNavneeth Basically, I separated the integrand into four quadrants and combined them symbolically before performing the integration. Doing so is very effective in this case, because of the substantial symbolic cancellations that occur. Thanks for asking this question, because it caused me to realize that I should have applied FullSimplify to the new integrand to achieve much greater numerical precision. Please see the new section of my answer. Commented Nov 26, 2019 at 22:09

I think you have a couple of typos in your formulas. Once these are corrected, the integral evaluates fine:

1. There should be no backtick at the end of the definition of Q[x], and
2. in the definition of F[x,y] I think you should replace (0.6329)^10^(-25) with (0.6329)*10^(-25).

Once these are fixed, the integral comes out as $$1.35767\times 10^{43}$$:

r = (1.082)*10^8;
h = (4.87)*10^(24);
Q[x_] = x - (0.206)*Sin[x] + (0.206)*(0.206)*Sin[x]*Cos[x];
R[x_] = ((57.91)*10^6)*(1 + 0.206*Cos[Q[x]]);
B[x_] = (0.206)*Sin[2*x] - Sin[x];
A[x_] = ((0.62)*10^(-25))*((4.854)*Cos[Q[x]] - B[Q[x]]*Sin[Q[x]]);
F[x_, y_] = R[x]*A[x] - r*Cos[y]*((((57.91)*10^6)*A[x]*(1 - (11.93)*10^6)*
Cos[Q[x]]/R[x]) - ((-3.008)*10^(-25)) - ((57.91)*10^6)*
B[x]*((-0.052)*10^(-31))*Sin[Q[x]]) -
r*Sin[y]*(((0.6329)*10^(-25))*
Sin[Q[x]] + ((56.67)*10^6)*((-0.05196)*10^(-31))*B[x]*
Cos[Q[x]]);
G[x_, y_] = (1 - (2*R[x]/r)*(Cos[y]*(((11.93)*10^6) + ((57.91)*10^6)*Cos[Q[x]])/
R[x]) + (R[x]/r)^2)^(-1.5);
u = ((0.01917)*10^(27))*((57.91)*10^6);
NIntegrate[(h/(((1.989)*10^(30)))*(r^3))*F[x, y]*G[x, y]*u,
{x, 0, 2 Pi}, {y, 0, 2 Pi}]

(*    1.35767*10^43    *)